Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-10-04 05:19:02 

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