Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 /**_________________________________________________________________
0002    class:   BSFitter.cc
0003    package: RecoVertex/BeamSpotProducer
0004 
0005 
0006 
0007  author: Francisco Yumiceva, Fermilab (yumiceva@fnal.gov)
0008 
0009 
0010 
0011 ________________________________________________________________**/
0012 
0013 #include "Minuit2/VariableMetricMinimizer.h"
0014 #include "Minuit2/FunctionMinimum.h"
0015 #include "Minuit2/MnPrint.h"
0016 #include "Minuit2/MnMigrad.h"
0017 #include "Minuit2/MnUserParameterState.h"
0018 //#include "CLHEP/config/CLHEP.h"
0019 
0020 // C++ standard
0021 #include <vector>
0022 #include <cmath>
0023 
0024 // CMS
0025 #include "RecoVertex/BeamSpotProducer/interface/BSFitter.h"
0026 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0027 #include "FWCore/Utilities/interface/Exception.h"
0028 #include "FWCore/Utilities/interface/isFinite.h"
0029 
0030 // ROOT
0031 #include "TMatrixD.h"
0032 #include "TMatrixDSym.h"
0033 #include "TDecompBK.h"
0034 #include "TH1.h"
0035 #include "TF1.h"
0036 #include "TMinuitMinimizer.h"
0037 
0038 using namespace ROOT::Minuit2;
0039 
0040 //_____________________________________________________________________
0041 BSFitter::BSFitter() {
0042   fbeamtype = reco::BeamSpot::Unknown;
0043 
0044   //In order to make fitting ROOT histograms thread safe
0045   // one must call this undocumented function
0046   TMinuitMinimizer::UseStaticMinuit(false);
0047 }
0048 
0049 //_____________________________________________________________________
0050 BSFitter::BSFitter(const std::vector<BSTrkParameters> &BSvector) {
0051   //In order to make fitting ROOT histograms thread safe
0052   // one must call this undocumented function
0053   TMinuitMinimizer::UseStaticMinuit(false);
0054 
0055   ffit_type = "default";
0056   ffit_variable = "default";
0057 
0058   fBSvector = BSvector;
0059 
0060   fsqrt2pi = sqrt(2. * TMath::Pi());
0061 
0062   fpar_name[0] = "z0        ";
0063   fpar_name[1] = "SigmaZ0   ";
0064   fpar_name[2] = "X0        ";
0065   fpar_name[3] = "Y0        ";
0066   fpar_name[4] = "dxdz      ";
0067   fpar_name[5] = "dydz      ";
0068   fpar_name[6] = "SigmaBeam ";
0069 
0070   //if (theGausszFcn == 0 ) {
0071   thePDF = new BSpdfsFcn();
0072 
0073   //}
0074   //if (theFitter == 0 ) {
0075 
0076   theFitter = new VariableMetricMinimizer();
0077 
0078   //}
0079 
0080   fapplyd0cut = false;
0081   fapplychi2cut = false;
0082   ftmprow = 0;
0083   ftmp.ResizeTo(4, 1);
0084   ftmp.Zero();
0085   fnthite = 0;
0086   fMaxZ = 50.;         //cm
0087   fconvergence = 0.5;  // stop fit when 50% of the input collection has been removed.
0088   fminNtrks = 100;
0089   finputBeamWidth = -1;  // no input
0090 
0091   h1z = new TH1F("h1z", "z distribution", 200, -fMaxZ, fMaxZ);
0092 }
0093 
0094 //______________________________________________________________________
0095 BSFitter::~BSFitter() {
0096   //delete fBSvector;
0097   delete thePDF;
0098   delete theFitter;
0099 }
0100 
0101 //______________________________________________________________________
0102 reco::BeamSpot BSFitter::Fit() { return this->Fit(nullptr); }
0103 
0104 //______________________________________________________________________
0105 reco::BeamSpot BSFitter::Fit(double *inipar = nullptr) {
0106   fbeamtype = reco::BeamSpot::Unknown;
0107   if (ffit_variable == "z") {
0108     if (ffit_type == "chi2") {
0109       return Fit_z_chi2(inipar);
0110 
0111     } else if (ffit_type == "likelihood") {
0112       return Fit_z_likelihood(inipar);
0113 
0114     } else if (ffit_type == "combined") {
0115       reco::BeamSpot tmp_beamspot = Fit_z_chi2(inipar);
0116       double tmp_par[2] = {tmp_beamspot.z0(), tmp_beamspot.sigmaZ()};
0117       return Fit_z_likelihood(tmp_par);
0118 
0119     } else {
0120       throw cms::Exception("LogicError")
0121           << "Error in BeamSpotProducer/BSFitter: "
0122           << "Illegal fit type, options are chi2,likelihood or combined(ie. first chi2 then likelihood)";
0123     }
0124   } else if (ffit_variable == "d") {
0125     if (ffit_type == "d0phi") {
0126       this->d0phi_Init();
0127       return Fit_d0phi();
0128 
0129     } else if (ffit_type == "likelihood") {
0130       return Fit_d_likelihood(inipar);
0131 
0132     } else if (ffit_type == "combined") {
0133       this->d0phi_Init();
0134       reco::BeamSpot tmp_beamspot = Fit_d0phi();
0135       double tmp_par[4] = {tmp_beamspot.x0(), tmp_beamspot.y0(), tmp_beamspot.dxdz(), tmp_beamspot.dydz()};
0136       return Fit_d_likelihood(tmp_par);
0137 
0138     } else {
0139       throw cms::Exception("LogicError") << "Error in BeamSpotProducer/BSFitter: "
0140                                          << "Illegal fit type, options are d0phi, likelihood or combined";
0141     }
0142   } else if (ffit_variable == "d*z" || ffit_variable == "default") {
0143     if (ffit_type == "likelihood" || ffit_type == "default") {
0144       reco::BeamSpot::CovarianceMatrix matrix;
0145       // we are now fitting Z inside d0phi fitter
0146       // first fit z distribution using a chi2 fit
0147       //reco::BeamSpot tmp_z = Fit_z_chi2(inipar);
0148       //for (int j = 2 ; j < 4 ; ++j) {
0149       //for(int k = j ; k < 4 ; ++k) {
0150       //    matrix(j,k) = tmp_z.covariance()(j,k);
0151       //}
0152       //}
0153 
0154       // use d0-phi algorithm to extract transverse position
0155       this->d0phi_Init();
0156       //reco::BeamSpot tmp_d0phi= Fit_d0phi(); // change to iterative procedure:
0157       this->Setd0Cut_d0phi(4.0);
0158       reco::BeamSpot tmp_d0phi = Fit_ited0phi();
0159 
0160       //for (int j = 0 ; j < 2 ; ++j) {
0161       //    for(int k = j ; k < 2 ; ++k) {
0162       //        matrix(j,k) = tmp_d0phi.covariance()(j,k);
0163       //}
0164       //}
0165       // slopes
0166       //for (int j = 4 ; j < 6 ; ++j) {
0167       // for(int k = j ; k < 6 ; ++k) {
0168       //  matrix(j,k) = tmp_d0phi.covariance()(j,k);
0169       //  }
0170       //}
0171 
0172       // put everything into one object
0173       reco::BeamSpot spot(reco::BeamSpot::Point(tmp_d0phi.x0(), tmp_d0phi.y0(), tmp_d0phi.z0()),
0174                           tmp_d0phi.sigmaZ(),
0175                           tmp_d0phi.dxdz(),
0176                           tmp_d0phi.dydz(),
0177                           0.,
0178                           tmp_d0phi.covariance(),
0179                           fbeamtype);
0180 
0181       //reco::BeamSpot tmp_z = Fit_z_chi2(inipar);
0182 
0183       //reco::BeamSpot tmp_d0phi = Fit_d0phi();
0184 
0185       // log-likelihood fit
0186       if (ffit_type == "likelihood") {
0187         double tmp_par[7] = {
0188             tmp_d0phi.x0(), tmp_d0phi.y0(), tmp_d0phi.z0(), tmp_d0phi.sigmaZ(), tmp_d0phi.dxdz(), tmp_d0phi.dydz(), 0.0};
0189 
0190         double tmp_error_par[7];
0191         for (int s = 0; s < 6; s++) {
0192           tmp_error_par[s] = pow(tmp_d0phi.covariance()(s, s), 0.5);
0193         }
0194         tmp_error_par[6] = 0.0;
0195 
0196         reco::BeamSpot tmp_lh = Fit_d_z_likelihood(tmp_par, tmp_error_par);
0197 
0198         if (edm::isNotFinite(ff_minimum)) {
0199           edm::LogWarning("BSFitter")
0200               << "BSFitter: Result is non physical. Log-Likelihood fit to extract beam width did not converge."
0201               << std::endl;
0202           tmp_lh.setType(reco::BeamSpot::Unknown);
0203           return tmp_lh;
0204         }
0205         return tmp_lh;
0206 
0207       } else {
0208         edm::LogInfo("BSFitter") << "default track-based fit does not extract beam width." << std::endl;
0209         return spot;
0210       }
0211 
0212     } else if (ffit_type == "resolution") {
0213       reco::BeamSpot tmp_z = Fit_z_chi2(inipar);
0214       this->d0phi_Init();
0215       reco::BeamSpot tmp_d0phi = Fit_d0phi();
0216 
0217       double tmp_par[7] = {
0218           tmp_d0phi.x0(), tmp_d0phi.y0(), tmp_z.z0(), tmp_z.sigmaZ(), tmp_d0phi.dxdz(), tmp_d0phi.dydz(), 0.0};
0219       double tmp_error_par[7];
0220       for (int s = 0; s < 6; s++) {
0221         tmp_error_par[s] = pow(tmp_par[s], 0.5);
0222       }
0223       tmp_error_par[6] = 0.0;
0224 
0225       reco::BeamSpot tmp_beam = Fit_d_z_likelihood(tmp_par, tmp_error_par);
0226 
0227       double tmp_par2[7] = {tmp_beam.x0(),
0228                             tmp_beam.y0(),
0229                             tmp_beam.z0(),
0230                             tmp_beam.sigmaZ(),
0231                             tmp_beam.dxdz(),
0232                             tmp_beam.dydz(),
0233                             tmp_beam.BeamWidthX()};
0234 
0235       reco::BeamSpot tmp_lh = Fit_dres_z_likelihood(tmp_par2);
0236 
0237       if (edm::isNotFinite(ff_minimum)) {
0238         edm::LogWarning("BSFitter") << "Result is non physical. Log-Likelihood fit did not converge." << std::endl;
0239         tmp_lh.setType(reco::BeamSpot::Unknown);
0240         return tmp_lh;
0241       }
0242       return tmp_lh;
0243 
0244     } else {
0245       throw cms::Exception("LogicError") << "Error in BeamSpotProducer/BSFitter: "
0246                                          << "Illegal fit type, options are likelihood or resolution";
0247     }
0248   } else {
0249     throw cms::Exception("LogicError") << "Error in BeamSpotProducer/BSFitter: "
0250                                        << "Illegal variable type, options are \"z\", \"d\", or \"d*z\"";
0251   }
0252 }
0253 
0254 //______________________________________________________________________
0255 reco::BeamSpot BSFitter::Fit_z_likelihood(double *inipar) {
0256   //std::cout << "Fit_z(double *) called" << std::endl;
0257   //std::cout << "inipar[0]= " << inipar[0] << std::endl;
0258   //std::cout << "inipar[1]= " << inipar[1] << std::endl;
0259 
0260   std::vector<double> par(2, 0);
0261   std::vector<double> err(2, 0);
0262 
0263   par.push_back(0.0);
0264   par.push_back(7.0);
0265   err.push_back(0.0001);
0266   err.push_back(0.0001);
0267   //par[0] = 0.0; err[0] = 0.0;
0268   //par[1] = 7.0; err[1] = 0.0;
0269 
0270   thePDF->SetPDFs("PDFGauss_z");
0271   thePDF->SetData(fBSvector);
0272   //std::cout << "data loaded"<< std::endl;
0273 
0274   //FunctionMinimum fmin = theFitter->Minimize(*theGausszFcn, par, err, 1, 500, 0.1);
0275   MnUserParameters upar;
0276   upar.Add("X0", 0., 0.);
0277   upar.Add("Y0", 0., 0.);
0278   upar.Add("Z0", inipar[0], 0.001);
0279   upar.Add("sigmaZ", inipar[1], 0.001);
0280 
0281   MnMigrad migrad(*thePDF, upar);
0282 
0283   FunctionMinimum fmin = migrad();
0284   ff_minimum = fmin.Fval();
0285   //std::cout << " eval= " << ff_minimum
0286   //          << "/n params[0]= " << fmin.Parameters().Vec()(0) << std::endl;
0287 
0288   /*
0289     TMinuit *gmMinuit = new TMinuit(2);
0290 
0291     //gmMinuit->SetFCN(z_fcn);
0292     gmMinuit->SetFCN(myFitz_fcn);
0293 
0294 
0295     int ierflg = 0;
0296     double step[2] = {0.001,0.001};
0297 
0298     for (int i = 0; i<2; i++) {
0299         gmMinuit->mnparm(i,fpar_name[i].c_str(),inipar[i],step[i],0,0,ierflg);
0300     }
0301     gmMinuit->Migrad();
0302     */
0303   reco::BeamSpot::CovarianceMatrix matrix;
0304 
0305   for (int j = 2; j < 4; ++j) {
0306     for (int k = j; k < 4; ++k) {
0307       matrix(j, k) = fmin.Error().Matrix()(j, k);
0308     }
0309   }
0310 
0311   return reco::BeamSpot(reco::BeamSpot::Point(0., 0., fmin.Parameters().Vec()(2)),
0312                         fmin.Parameters().Vec()(3),
0313                         0.,
0314                         0.,
0315                         0.,
0316                         matrix,
0317                         fbeamtype);
0318 }
0319 
0320 //______________________________________________________________________
0321 reco::BeamSpot BSFitter::Fit_z_chi2(double *inipar) {
0322   // N.B. this fit is not performed anymore but now
0323   // Z is fitted in the same track set used in the d0-phi fit after
0324   // each iteration
0325 
0326   //std::cout << "Fit_z_chi2() called" << std::endl;
0327   // FIXME: include whole tracker z length for the time being
0328   // ==> add protection and z0 cut
0329   h1z = new TH1F("h1z", "z distribution", 200, -fMaxZ, fMaxZ);
0330 
0331   std::vector<BSTrkParameters>::const_iterator iparam = fBSvector.begin();
0332 
0333   // HERE check size of track vector
0334 
0335   for (iparam = fBSvector.begin(); iparam != fBSvector.end(); ++iparam) {
0336     h1z->Fill(iparam->z0());
0337     //std::cout<<"z0="<<iparam->z0()<<"; sigZ0="<<iparam->sigz0()<<std::endl;
0338   }
0339 
0340   //Use our own copy for thread safety
0341   // also do not add to global list of functions
0342   TF1 fgaus("fgaus", "gaus", 0., 1., TF1::EAddToList::kNo);
0343   h1z->Fit(&fgaus, "QLMN0 SERIAL");
0344   //std::cout << "fitted "<< std::endl;
0345 
0346   //std::cout << "got function" << std::endl;
0347   double fpar[2] = {fgaus.GetParameter(1), fgaus.GetParameter(2)};
0348   //std::cout<<"Debug fpar[2] = (" <<fpar[0]<<","<<fpar[1]<<")"<<std::endl;
0349   reco::BeamSpot::CovarianceMatrix matrix;
0350   // add matrix values.
0351   matrix(2, 2) = fgaus.GetParError(1) * fgaus.GetParError(1);
0352   matrix(3, 3) = fgaus.GetParError(2) * fgaus.GetParError(2);
0353 
0354   //delete h1z;
0355 
0356   return reco::BeamSpot(reco::BeamSpot::Point(0., 0., fpar[0]), fpar[1], 0., 0., 0., matrix, fbeamtype);
0357 }
0358 
0359 //______________________________________________________________________
0360 reco::BeamSpot BSFitter::Fit_ited0phi() {
0361   this->d0phi_Init();
0362   edm::LogInfo("BSFitter") << "number of total input tracks: " << fBSvector.size() << std::endl;
0363 
0364   reco::BeamSpot theanswer;
0365 
0366   if ((int)fBSvector.size() <= fminNtrks) {
0367     edm::LogWarning("BSFitter") << "need at least " << fminNtrks << " tracks to run beamline fitter." << std::endl;
0368     fbeamtype = reco::BeamSpot::Fake;
0369     theanswer.setType(fbeamtype);
0370     return theanswer;
0371   }
0372 
0373   theanswer = Fit_d0phi();  //get initial ftmp and ftmprow
0374   if (goodfit)
0375     fnthite++;
0376   //std::cout << "Initial tempanswer (iteration 0): " << theanswer << std::endl;
0377 
0378   reco::BeamSpot preanswer = theanswer;
0379 
0380   while (goodfit && ftmprow > fconvergence * fBSvector.size() && ftmprow > fminNtrks) {
0381     theanswer = Fit_d0phi();
0382     fd0cut /= 1.5;
0383     fchi2cut /= 1.5;
0384     if (goodfit && ftmprow > fconvergence * fBSvector.size() && ftmprow > fminNtrks) {
0385       preanswer = theanswer;
0386       //std::cout << "Iteration " << fnthite << ": " << preanswer << std::endl;
0387       fnthite++;
0388     }
0389   }
0390   // FIXME: return fit results from previous iteration for both bad fit and for >50% tracks thrown away
0391   //std::cout << "The last iteration, theanswer: " << theanswer << std::endl;
0392   theanswer = preanswer;
0393   //std::cout << "Use previous results from iteration #" << ( fnthite > 0 ? fnthite-1 : 0 ) << std::endl;
0394   //if ( fnthite > 1 ) std::cout << theanswer << std::endl;
0395 
0396   edm::LogInfo("BSFitter") << "Total number of successful iterations = " << (goodfit ? (fnthite + 1) : fnthite)
0397                            << std::endl;
0398   if (goodfit) {
0399     fbeamtype = reco::BeamSpot::Tracker;
0400     theanswer.setType(fbeamtype);
0401   } else {
0402     edm::LogWarning("BSFitter") << "Fit doesn't converge!!!" << std::endl;
0403     fbeamtype = reco::BeamSpot::Unknown;
0404     theanswer.setType(fbeamtype);
0405   }
0406   return theanswer;
0407 }
0408 
0409 //______________________________________________________________________
0410 reco::BeamSpot BSFitter::Fit_d0phi() {
0411   //LogDebug ("BSFitter") << " we will use " << fBSvector.size() << " tracks.";
0412   if (fnthite > 0)
0413     edm::LogInfo("BSFitter") << " number of tracks used: " << ftmprow << std::endl;
0414   //std::cout << " ftmp = matrix("<<ftmp.GetNrows()<<","<<ftmp.GetNcols()<<")"<<std::endl;
0415   //std::cout << " ftmp(0,0)="<<ftmp(0,0)<<std::endl;
0416   //std::cout << " ftmp(1,0)="<<ftmp(1,0)<<std::endl;
0417   //std::cout << " ftmp(2,0)="<<ftmp(2,0)<<std::endl;
0418   //std::cout << " ftmp(3,0)="<<ftmp(3,0)<<std::endl;
0419 
0420   h1z->Reset();
0421 
0422   TMatrixD x_result(4, 1);
0423   TMatrixDSym V_result(4);
0424 
0425   TMatrixDSym Vint(4);
0426   TMatrixD b(4, 1);
0427 
0428   //Double_t weightsum = 0;
0429 
0430   Vint.Zero();
0431   b.Zero();
0432 
0433   TMatrixD g(4, 1);
0434   TMatrixDSym temp(4);
0435 
0436   std::vector<BSTrkParameters>::iterator iparam = fBSvector.begin();
0437   ftmprow = 0;
0438 
0439   //edm::LogInfo ("BSFitter") << " test";
0440 
0441   //std::cout << "BSFitter: fit" << std::endl;
0442 
0443   for (iparam = fBSvector.begin(); iparam != fBSvector.end(); ++iparam) {
0444     //if(i->weight2 == 0) continue;
0445 
0446     //if (ftmprow==0) {
0447     //std::cout << "d0=" << iparam->d0() << " sigd0=" << iparam->sigd0()
0448     //<< " phi0="<< iparam->phi0() << " z0=" << iparam->z0() << std::endl;
0449     //std::cout << "d0phi_d0=" << iparam->d0phi_d0() << " d0phi_chi2="<<iparam->d0phi_chi2() << std::endl;
0450     //}
0451     g(0, 0) = sin(iparam->phi0());
0452     g(1, 0) = -cos(iparam->phi0());
0453     g(2, 0) = iparam->z0() * g(0, 0);
0454     g(3, 0) = iparam->z0() * g(1, 0);
0455 
0456     // average transverse beam width
0457     double sigmabeam2 = 0.006 * 0.006;
0458     if (finputBeamWidth > 0)
0459       sigmabeam2 = finputBeamWidth * finputBeamWidth;
0460     else {
0461       //edm::LogWarning("BSFitter") << "using in fit beam width = " << sqrt(sigmabeam2) << std::endl;
0462     }
0463 
0464     //double sigma2 = sigmabeam2 +  (iparam->sigd0())* (iparam->sigd0()) / iparam->weight2;
0465     // this should be 2*sigmabeam2?
0466     double sigma2 = sigmabeam2 + (iparam->sigd0()) * (iparam->sigd0());
0467 
0468     TMatrixD ftmptrans(1, 4);
0469     ftmptrans = ftmptrans.Transpose(ftmp);
0470     TMatrixD dcor = ftmptrans * g;
0471     double chi2tmp = (iparam->d0() - dcor(0, 0)) * (iparam->d0() - dcor(0, 0)) / sigma2;
0472     (*iparam) = BSTrkParameters(
0473         iparam->z0(), iparam->sigz0(), iparam->d0(), iparam->sigd0(), iparam->phi0(), iparam->pt(), dcor(0, 0), chi2tmp);
0474 
0475     bool pass = true;
0476     if (fapplyd0cut && fnthite > 0) {
0477       if (std::abs(iparam->d0() - dcor(0, 0)) > fd0cut)
0478         pass = false;
0479     }
0480     if (fapplychi2cut && fnthite > 0) {
0481       if (chi2tmp > fchi2cut)
0482         pass = false;
0483     }
0484 
0485     if (pass) {
0486       temp.Zero();
0487       for (int j = 0; j < 4; ++j) {
0488         for (int k = j; k < 4; ++k) {
0489           temp(j, k) += g(j, 0) * g(k, 0);
0490         }
0491       }
0492 
0493       Vint += (temp * (1 / sigma2));
0494       b += (iparam->d0() / sigma2 * g);
0495       //weightsum += sqrt(i->weight2);
0496       ftmprow++;
0497       h1z->Fill(iparam->z0());
0498     }
0499   }
0500   Double_t determinant;
0501   TDecompBK bk(Vint);
0502   bk.SetTol(1e-11);  //FIXME: find a better way to solve x_result
0503   if (!bk.Decompose()) {
0504     goodfit = false;
0505     edm::LogWarning("BSFitter") << "Decomposition failed, matrix singular ?" << std::endl
0506                                 << "condition number = " << bk.Condition() << std::endl;
0507   } else {
0508     V_result = Vint.InvertFast(&determinant);
0509     x_result = V_result * b;
0510   }
0511   //        for(int j = 0 ; j < 4 ; ++j) {
0512   //      for(int k = 0 ; k < 4 ; ++k) {
0513   //        std::cout<<"V_result("<<j<<","<<k<<")="<<V_result(j,k)<<std::endl;
0514   //      }
0515   //        }
0516   //         for (int j=0;j<4;++j){
0517   //      std::cout<<"x_result("<<j<<",0)="<<x_result(j,0)<<std::endl;
0518   //    }
0519   //LogDebug ("BSFitter") << " d0-phi fit done.";
0520   //std::cout<< " d0-phi fit done." << std::endl;
0521 
0522   //Use our own copy for thread safety
0523   // also do not add to global list of functions
0524   TF1 fgaus("fgaus", "gaus", 0., 1., TF1::EAddToList::kNo);
0525   //returns 0 if OK
0526   //auto status = h1z->Fit(&fgaus,"QLM0","",h1z->GetMean() -2.*h1z->GetRMS(),h1z->GetMean() +2.*h1z->GetRMS());
0527   auto status =
0528       h1z->Fit(&fgaus, "QLN0 SERIAL", "", h1z->GetMean() - 2. * h1z->GetRMS(), h1z->GetMean() + 2. * h1z->GetRMS());
0529 
0530   //std::cout << "fitted "<< std::endl;
0531 
0532   //std::cout << "got function" << std::endl;
0533   if (status) {
0534     //edm::LogError("NoBeamSpotFit")<<"gaussian fit failed. no BS d0 fit";
0535 
0536     goodfit = false;
0537     return reco::BeamSpot();
0538   }
0539   double fpar[2] = {fgaus.GetParameter(1), fgaus.GetParameter(2)};
0540 
0541   reco::BeamSpot::CovarianceMatrix matrix;
0542   // first two parameters
0543   for (int j = 0; j < 2; ++j) {
0544     for (int k = j; k < 2; ++k) {
0545       matrix(j, k) = V_result(j, k);
0546     }
0547   }
0548   // slope parameters
0549   for (int j = 4; j < 6; ++j) {
0550     for (int k = j; k < 6; ++k) {
0551       matrix(j, k) = V_result(j - 2, k - 2);
0552     }
0553   }
0554 
0555   // Z0 and sigmaZ
0556   matrix(2, 2) = fgaus.GetParError(1) * fgaus.GetParError(1);
0557   matrix(3, 3) = fgaus.GetParError(2) * fgaus.GetParError(2);
0558 
0559   ftmp = x_result;
0560 
0561   // x0 and y0 are *not* x,y at z=0, but actually at z=0
0562   // to correct for this, we need to translate them to z=z0
0563   // using the measured slopes
0564   //
0565   double x0tmp = x_result(0, 0);
0566   double y0tmp = x_result(1, 0);
0567 
0568   x0tmp += x_result(2, 0) * fpar[0];
0569   y0tmp += x_result(3, 0) * fpar[0];
0570 
0571   return reco::BeamSpot(
0572       reco::BeamSpot::Point(x0tmp, y0tmp, fpar[0]), fpar[1], x_result(2, 0), x_result(3, 0), 0., matrix, fbeamtype);
0573 }
0574 
0575 //______________________________________________________________________
0576 void BSFitter::Setd0Cut_d0phi(double d0cut) {
0577   fapplyd0cut = true;
0578 
0579   //fBSforCuts = BSfitted;
0580   fd0cut = d0cut;
0581 }
0582 
0583 //______________________________________________________________________
0584 void BSFitter::SetChi2Cut_d0phi(double chi2cut) {
0585   fapplychi2cut = true;
0586 
0587   //fBSforCuts = BSfitted;
0588   fchi2cut = chi2cut;
0589 }
0590 
0591 //______________________________________________________________________
0592 reco::BeamSpot BSFitter::Fit_d_likelihood(double *inipar) {
0593   thePDF->SetPDFs("PDFGauss_d");
0594   thePDF->SetData(fBSvector);
0595 
0596   MnUserParameters upar;
0597   upar.Add("X0", inipar[0], 0.001);
0598   upar.Add("Y0", inipar[1], 0.001);
0599   upar.Add("Z0", 0., 0.001);
0600   upar.Add("sigmaZ", 0., 0.001);
0601   upar.Add("dxdz", inipar[2], 0.001);
0602   upar.Add("dydz", inipar[3], 0.001);
0603 
0604   MnMigrad migrad(*thePDF, upar);
0605 
0606   FunctionMinimum fmin = migrad();
0607   ff_minimum = fmin.Fval();
0608 
0609   reco::BeamSpot::CovarianceMatrix matrix;
0610   for (int j = 0; j < 6; ++j) {
0611     for (int k = j; k < 6; ++k) {
0612       matrix(j, k) = fmin.Error().Matrix()(j, k);
0613     }
0614   }
0615 
0616   return reco::BeamSpot(reco::BeamSpot::Point(fmin.Parameters().Vec()(0), fmin.Parameters().Vec()(1), 0.),
0617                         0.,
0618                         fmin.Parameters().Vec()(4),
0619                         fmin.Parameters().Vec()(5),
0620                         0.,
0621                         matrix,
0622                         fbeamtype);
0623 }
0624 //______________________________________________________________________
0625 double BSFitter::scanPDF(double *init_pars, int &tracksfixed, int option) {
0626   if (option == 1)
0627     init_pars[6] = 0.0005;  //starting value for any given configuration
0628 
0629   //local vairables with initial values
0630   double fsqrt2pi = 0.0;
0631   double d_sig = 0.0;
0632   double d_dprime = 0.0;
0633   double d_result = 0.0;
0634   double z_sig = 0.0;
0635   double z_result = 0.0;
0636   double function = 0.0;
0637   double tot_pdf = 0.0;
0638   double last_minvalue = 1.0e+10;
0639   double init_bw = -99.99;
0640   int iters = 0;
0641 
0642   //used to remove tracks if far away from bs by this
0643   double DeltadCut = 0.1000;
0644   if (init_pars[6] < 0.0200) {
0645     DeltadCut = 0.0900;
0646   }  //worked for high 2.36TeV
0647   if (init_pars[6] < 0.0100) {
0648     DeltadCut = 0.0700;
0649   }  //just a guesss for 7 TeV but one should scan for actual values
0650 
0651   std::vector<BSTrkParameters>::const_iterator iparam = fBSvector.begin();
0652 
0653   if (option == 1)
0654     iters = 500;
0655   if (option == 2)
0656     iters = 1;
0657 
0658   for (int p = 0; p < iters; p++) {
0659     if (iters == 500)
0660       init_pars[6] += 0.0002;
0661     tracksfixed = 0;
0662 
0663     for (iparam = fBSvector.begin(); iparam != fBSvector.end(); ++iparam) {
0664       fsqrt2pi = sqrt(2. * TMath::Pi());
0665       d_sig = sqrt(init_pars[6] * init_pars[6] + (iparam->sigd0()) * (iparam->sigd0()));
0666       d_dprime = iparam->d0() - (((init_pars[0] + iparam->z0() * (init_pars[4])) * sin(iparam->phi0())) -
0667                                  ((init_pars[1] + iparam->z0() * (init_pars[5])) * cos(iparam->phi0())));
0668 
0669       //***Remove tracks before the fit which gives low pdf values to blow up the pdf
0670       if (std::abs(d_dprime) < DeltadCut && option == 2) {
0671         fBSvectorBW.push_back(*iparam);
0672       }
0673 
0674       d_result = (exp(-(d_dprime * d_dprime) / (2.0 * d_sig * d_sig))) / (d_sig * fsqrt2pi);
0675       z_sig = sqrt(iparam->sigz0() * iparam->sigz0() + init_pars[3] * init_pars[3]);
0676       z_result = (exp(-((iparam->z0() - init_pars[2]) * (iparam->z0() - init_pars[2])) / (2.0 * z_sig * z_sig))) /
0677                  (z_sig * fsqrt2pi);
0678       tot_pdf = z_result * d_result;
0679 
0680       //for those trcks which gives problems due to very tiny pdf_d values.
0681       //Update: This protection will NOT be used with the dprime cut above but still kept here to get
0682       // the intial value of beam width reasonably
0683       //A warning will appear if there were any tracks with < 10^-5 for pdf_d so that (d-dprime) cut can be lowered
0684       if (d_result < 1.0e-05) {
0685         tot_pdf = z_result * 1.0e-05;
0686         //if(option==2)std::cout<<"last Iter  d-d'   =  "<<(std::abs(d_dprime))<<std::endl;
0687         tracksfixed++;
0688       }
0689 
0690       function = function + log(tot_pdf);
0691 
0692     }  //loop over tracks
0693 
0694     function = -2.0 * function;
0695     if (function < last_minvalue) {
0696       init_bw = init_pars[6];
0697       last_minvalue = function;
0698     }
0699     function = 0.0;
0700   }  //loop over beam width
0701 
0702   if (init_bw > 0) {
0703     init_bw = init_bw + (0.20 * init_bw);  //start with 20 % more
0704 
0705   } else {
0706     if (option == 1) {
0707       edm::LogWarning("BSFitter")
0708           << "scanPDF:====>>>> WARNING***: The initial guess value of Beam width is negative!!!!!!" << std::endl
0709           << "scanPDF:====>>>> Assigning beam width a starting value of " << init_bw << "  cm" << std::endl;
0710       init_bw = 0.0200;
0711     }
0712   }
0713 
0714   return init_bw;
0715 }
0716 
0717 //________________________________________________________________________________
0718 reco::BeamSpot BSFitter::Fit_d_z_likelihood(double *inipar, double *error_par) {
0719   int tracksFailed = 0;
0720 
0721   //estimate first guess of beam width and tame 20% extra of it to start
0722   inipar[6] = scanPDF(inipar, tracksFailed, 1);
0723   error_par[6] = (inipar[6]) * 0.20;
0724 
0725   //Here remove the tracks which give low pdf and fill into a new vector
0726   //std::cout<<"Size of Old vector = "<<(fBSvector.size())<<std::endl;
0727   /* double junk= */ scanPDF(inipar, tracksFailed, 2);
0728   //std::cout<<"Size of New vector = "<<(fBSvectorBW.size())<<std::endl;
0729 
0730   //Refill the fBSVector again with new sets of tracks
0731   fBSvector.clear();
0732   std::vector<BSTrkParameters>::const_iterator iparamBW = fBSvectorBW.begin();
0733   for (iparamBW = fBSvectorBW.begin(); iparamBW != fBSvectorBW.end(); ++iparamBW) {
0734     fBSvector.push_back(*iparamBW);
0735   }
0736 
0737   thePDF->SetPDFs("PDFGauss_d*PDFGauss_z");
0738   thePDF->SetData(fBSvector);
0739   MnUserParameters upar;
0740 
0741   upar.Add("X0", inipar[0], error_par[0]);
0742   upar.Add("Y0", inipar[1], error_par[1]);
0743   upar.Add("Z0", inipar[2], error_par[2]);
0744   upar.Add("sigmaZ", inipar[3], error_par[3]);
0745   upar.Add("dxdz", inipar[4], error_par[4]);
0746   upar.Add("dydz", inipar[5], error_par[5]);
0747   upar.Add("BeamWidthX", inipar[6], error_par[6]);
0748 
0749   MnMigrad migrad(*thePDF, upar);
0750 
0751   FunctionMinimum fmin = migrad();
0752 
0753   // std::cout<<"-----how the fit evoves------"<<std::endl;
0754   // std::cout<<fmin<<std::endl;
0755 
0756   ff_minimum = fmin.Fval();
0757 
0758   bool ff_nfcn = fmin.HasReachedCallLimit();
0759   bool ff_cov = fmin.HasCovariance();
0760   bool testing = fmin.IsValid();
0761 
0762   //Print WARNINGS if minimum did not converged
0763   if (!testing) {
0764     edm::LogWarning("BSFitter") << "===========>>>>>** WARNING: MINUIT DID NOT CONVERGES PROPERLY !!!!!!" << std::endl;
0765     if (ff_nfcn)
0766       edm::LogWarning("BSFitter") << "===========>>>>>** WARNING: No. of Calls Exhausted" << std::endl;
0767     if (!ff_cov)
0768       edm::LogWarning("BSFitter") << "===========>>>>>** WARNING: Covariance did not found" << std::endl;
0769   }
0770 
0771   edm::LogInfo("BSFitter") << "The Total # Tracks used for beam width fit = " << (fBSvectorBW.size()) << std::endl;
0772 
0773   //Checks after fit is performed
0774   double lastIter_pars[7];
0775 
0776   for (int ip = 0; ip < 7; ip++) {
0777     lastIter_pars[ip] = fmin.Parameters().Vec()(ip);
0778   }
0779 
0780   tracksFailed = 0;
0781   /* double lastIter_scan= */ scanPDF(lastIter_pars, tracksFailed, 2);
0782 
0783   edm::LogWarning("BSFitter") << "WARNING: # of tracks which have very low pdf value (pdf_d < 1.0e-05) are  = "
0784                               << tracksFailed << std::endl;
0785 
0786   //std::cout << " eval= " << ff_minimum
0787   //                << "/n params[0]= " << fmin.Parameters().Vec()(0) << std::endl;
0788 
0789   reco::BeamSpot::CovarianceMatrix matrix;
0790 
0791   for (int j = 0; j < 7; ++j) {
0792     for (int k = j; k < 7; ++k) {
0793       matrix(j, k) = fmin.Error().Matrix()(j, k);
0794     }
0795   }
0796 
0797   return reco::BeamSpot(
0798       reco::BeamSpot::Point(fmin.Parameters().Vec()(0), fmin.Parameters().Vec()(1), fmin.Parameters().Vec()(2)),
0799       fmin.Parameters().Vec()(3),
0800       fmin.Parameters().Vec()(4),
0801       fmin.Parameters().Vec()(5),
0802       fmin.Parameters().Vec()(6),
0803 
0804       matrix,
0805       fbeamtype);
0806 }
0807 
0808 //______________________________________________________________________
0809 reco::BeamSpot BSFitter::Fit_dres_z_likelihood(double *inipar) {
0810   thePDF->SetPDFs("PDFGauss_d_resolution*PDFGauss_z");
0811   thePDF->SetData(fBSvector);
0812 
0813   MnUserParameters upar;
0814   upar.Add("X0", inipar[0], 0.001);
0815   upar.Add("Y0", inipar[1], 0.001);
0816   upar.Add("Z0", inipar[2], 0.001);
0817   upar.Add("sigmaZ", inipar[3], 0.001);
0818   upar.Add("dxdz", inipar[4], 0.001);
0819   upar.Add("dydz", inipar[5], 0.001);
0820   upar.Add("BeamWidthX", inipar[6], 0.0001);
0821   upar.Add("c0", 0.0010, 0.0001);
0822   upar.Add("c1", 0.0090, 0.0001);
0823 
0824   // fix beam width
0825   upar.Fix("BeamWidthX");
0826   // number of parameters in fit are 9-1 = 8
0827 
0828   MnMigrad migrad(*thePDF, upar);
0829 
0830   FunctionMinimum fmin = migrad();
0831   ff_minimum = fmin.Fval();
0832 
0833   reco::BeamSpot::CovarianceMatrix matrix;
0834 
0835   for (int j = 0; j < 6; ++j) {
0836     for (int k = j; k < 6; ++k) {
0837       matrix(j, k) = fmin.Error().Matrix()(j, k);
0838     }
0839   }
0840 
0841   //std::cout << " fill resolution values" << std::endl;
0842   //std::cout << " matrix size= " << fmin.Error().Matrix().size() << std::endl;
0843   //std::cout << " vec(6)="<< fmin.Parameters().Vec()(6) << std::endl;
0844   //std::cout << " vec(7)="<< fmin.Parameters().Vec()(7) << std::endl;
0845 
0846   fresolution_c0 = fmin.Parameters().Vec()(6);
0847   fresolution_c1 = fmin.Parameters().Vec()(7);
0848   fres_c0_err = sqrt(fmin.Error().Matrix()(6, 6));
0849   fres_c1_err = sqrt(fmin.Error().Matrix()(7, 7));
0850 
0851   for (int j = 6; j < 8; ++j) {
0852     for (int k = 6; k < 8; ++k) {
0853       fres_matrix(j - 6, k - 6) = fmin.Error().Matrix()(j, k);
0854     }
0855   }
0856 
0857   return reco::BeamSpot(
0858       reco::BeamSpot::Point(fmin.Parameters().Vec()(0), fmin.Parameters().Vec()(1), fmin.Parameters().Vec()(2)),
0859       fmin.Parameters().Vec()(3),
0860       fmin.Parameters().Vec()(4),
0861       fmin.Parameters().Vec()(5),
0862       inipar[6],
0863       matrix,
0864       fbeamtype);
0865 }
[ 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!