Reference Manual

Go to the documentation of this file.
00001 
00015 #include "Minuit2/VariableMetricMinimizer.h"
00016 #include "Minuit2/FunctionMinimum.h"
00017 #include "Minuit2/MnPrint.h"
00018 #include "Minuit2/MnMigrad.h"
00019 #include "Minuit2/MnUserParameterState.h"
00020 //#include "CLHEP/config/CLHEP.h"
00021 
00022 // C++ standard
00023 #include <vector>
00024 #include <cmath>
00025 
00026 // CMS
00027 #include "RecoVertex/BeamSpotProducer/interface/BSFitter.h"
00028 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00029 #include "FWCore/Utilities/interface/Exception.h"
00030 
00031 // ROOT
00032 #include "TMatrixD.h"
00033 #include "TMatrixDSym.h"
00034 #include "TDecompBK.h"
00035 #include "TH1.h"
00036 #include "TF1.h"
00037 
00038 using namespace ROOT::Minuit2;
00039 
00040 
00041 //_____________________________________________________________________
00042 BSFitter::BSFitter() {
00043         fbeamtype = reco::BeamSpot::Unknown;
00044 }
00045 
00046 //_____________________________________________________________________
00047 BSFitter::BSFitter( std:: vector< BSTrkParameters > BSvector ) {
00048 
00049         ffit_type = "default";
00050         ffit_variable = "default";
00051         
00052         fBSvector = BSvector;
00053 
00054         fsqrt2pi = sqrt(2.* TMath::Pi());
00055         
00056         fpar_name[0] = "z0        ";
00057         fpar_name[1] = "SigmaZ0   ";
00058         fpar_name[2] = "X0        ";
00059         fpar_name[3] = "Y0        ";
00060         fpar_name[4] = "dxdz      ";
00061         fpar_name[5] = "dydz      ";
00062         fpar_name[6] = "SigmaBeam ";               
00063 
00064         //if (theGausszFcn == 0 ) {
00065         thePDF = new BSpdfsFcn();
00066                 
00067 
00068 //}
00069                 //if (theFitter == 0 ) {
00070                 
00071         theFitter    = new VariableMetricMinimizer();
00072                 
00073                 //}
00074     
00075         fapplyd0cut = false;
00076         fapplychi2cut = false;
00077         ftmprow = 0;
00078         ftmp.ResizeTo(4,1);
00079         ftmp.Zero();
00080         fnthite=0;
00081         fMaxZ = 50.; //cm
00082         fconvergence = 0.5; // stop fit when 50% of the input collection has been removed.
00083         fminNtrks = 100;
00084         finputBeamWidth = -1; // no input
00085 
00086     h1z = new TH1F("h1z","z distribution",200,-fMaxZ, fMaxZ);
00087         
00088 }
00089 
00090 //______________________________________________________________________
00091 BSFitter::~BSFitter()
00092 {
00093         //delete fBSvector;
00094         delete thePDF;
00095         delete theFitter;
00096 }
00097 
00098 
00099 //______________________________________________________________________
00100 reco::BeamSpot BSFitter::Fit() {
00101 
00102         return this->Fit(0);
00103         
00104 }
00105 
00106 //______________________________________________________________________
00107 reco::BeamSpot BSFitter::Fit(double *inipar = 0) {
00108         fbeamtype = reco::BeamSpot::Unknown;
00109         if ( ffit_variable == "z" ) {
00110 
00111                 if ( ffit_type == "chi2" ) {
00112 
00113                         return Fit_z_chi2(inipar);
00114                         
00115                 } else if ( ffit_type == "likelihood" ) {
00116 
00117                         return Fit_z_likelihood(inipar);
00118                         
00119                 } else if ( ffit_type == "combined" ) {
00120 
00121                         reco::BeamSpot tmp_beamspot = Fit_z_chi2(inipar);
00122                         double tmp_par[2] = {tmp_beamspot.z0(), tmp_beamspot.sigmaZ()};
00123                         return Fit_z_likelihood(tmp_par);
00124                         
00125                 } else {
00126 
00127                         throw cms::Exception("LogicError")
00128                         << "Error in BeamSpotProducer/BSFitter: "
00129                         << "Illegal fit type, options are chi2,likelihood or combined(ie. first chi2 then likelihood)";
00130                         
00131                 }               
00132         } else if ( ffit_variable == "d" ) {
00133 
00134                 if ( ffit_type == "d0phi" ) {
00135                         this->d0phi_Init();
00136                         return Fit_d0phi();
00137                         
00138                 } else if ( ffit_type == "likelihood" ) {
00139 
00140                         return Fit_d_likelihood(inipar);
00141                         
00142                 } else if ( ffit_type == "combined" ) {
00143 
00144                         this->d0phi_Init();
00145                         reco::BeamSpot tmp_beamspot = Fit_d0phi();
00146                         double tmp_par[4] = {tmp_beamspot.x0(), tmp_beamspot.y0(), tmp_beamspot.dxdz(), tmp_beamspot.dydz()};
00147                         return Fit_d_likelihood(tmp_par);
00148                         
00149                 } else {
00150                         throw cms::Exception("LogicError")
00151                                 << "Error in BeamSpotProducer/BSFitter: "
00152                                 << "Illegal fit type, options are d0phi, likelihood or combined";
00153                 }
00154         } else if ( ffit_variable == "d*z" || ffit_variable == "default" ) {
00155 
00156                 if ( ffit_type == "likelihood" || ffit_type == "default" ) {
00157 
00158                         reco::BeamSpot::CovarianceMatrix matrix;
00159             // we are now fitting Z inside d0phi fitter 
00160                         // first fit z distribution using a chi2 fit
00161                         //reco::BeamSpot tmp_z = Fit_z_chi2(inipar);
00162                         //for (int j = 2 ; j < 4 ; ++j) {
00163             //for(int k = j ; k < 4 ; ++k) {
00164             //  matrix(j,k) = tmp_z.covariance()(j,k);
00165             //}
00166                         //}
00167                 
00168                         // use d0-phi algorithm to extract transverse position
00169                         this->d0phi_Init();
00170                         //reco::BeamSpot tmp_d0phi= Fit_d0phi(); // change to iterative procedure:
00171                         this->Setd0Cut_d0phi(4.0);
00172                         reco::BeamSpot tmp_d0phi= Fit_ited0phi();
00173                         
00174                         //for (int j = 0 ; j < 2 ; ++j) {
00175                         //      for(int k = j ; k < 2 ; ++k) {
00176                         //              matrix(j,k) = tmp_d0phi.covariance()(j,k);
00177             //}
00178                         //}
00179                         // slopes
00180                         //for (int j = 4 ; j < 6 ; ++j) {
00181             // for(int k = j ; k < 6 ; ++k) {
00182             //  matrix(j,k) = tmp_d0phi.covariance()(j,k);
00183                         //  }
00184                         //}
00185 
00186                 
00187                         // put everything into one object
00188                         reco::BeamSpot spot(reco::BeamSpot::Point(tmp_d0phi.x0(), tmp_d0phi.y0(), tmp_d0phi.z0()),
00189                                                                 tmp_d0phi.sigmaZ(),
00190                                                                 tmp_d0phi.dxdz(),
00191                                                                 tmp_d0phi.dydz(),
00192                                                                 0.,
00193                                                                 tmp_d0phi.covariance(),
00194                                                                 fbeamtype );
00195 
00196 
00197                         
00198                         //reco::BeamSpot tmp_z = Fit_z_chi2(inipar);
00199                         
00200                         //reco::BeamSpot tmp_d0phi = Fit_d0phi();
00201 
00202             // log-likelihood fit          
00203                         if (ffit_type == "likelihood") {
00204                 double tmp_par[7] = {tmp_d0phi.x0(), tmp_d0phi.y0(), tmp_d0phi.z0(),
00205                                      tmp_d0phi.sigmaZ(), tmp_d0phi.dxdz(), tmp_d0phi.dydz(),0.0};
00206                 
00207                 double tmp_error_par[7];
00208                 for(int s=0;s<6;s++){ tmp_error_par[s] = pow( tmp_d0phi.covariance()(s,s),0.5);}
00209                 tmp_error_par[6]=0.0;
00210                 
00211                 reco::BeamSpot tmp_lh = Fit_d_z_likelihood(tmp_par,tmp_error_par);
00212                 
00213                 if ( isnan(ff_minimum) || std::isinf(ff_minimum) ) {
00214                     edm::LogWarning("BSFitter") << "BSFitter: Result is non physical. Log-Likelihood fit to extract beam width did not converge." << std::endl;
00215                     tmp_lh.setType(reco::BeamSpot::Unknown);
00216                     return tmp_lh;                    
00217                 }
00218                 return tmp_lh;
00219                 
00220                         } else {
00221             
00222                 edm::LogInfo("BSFitter") << "default track-based fit does not extract beam width." << std::endl;
00223                                 return spot;
00224             }
00225                         
00226                         
00227                 } else if ( ffit_type == "resolution" ) {
00228 
00229                         reco::BeamSpot tmp_z = Fit_z_chi2(inipar);
00230                         this->d0phi_Init();                     
00231                         reco::BeamSpot tmp_d0phi = Fit_d0phi();
00232                         
00233                         double tmp_par[7] = {tmp_d0phi.x0(), tmp_d0phi.y0(), tmp_z.z0(),
00234                                                                  tmp_z.sigmaZ(), tmp_d0phi.dxdz(), tmp_d0phi.dydz(),0.0};
00235             double tmp_error_par[7];
00236             for(int s=0;s<6;s++){ tmp_error_par[s] = pow(tmp_par[s],0.5);}
00237             tmp_error_par[6]=0.0;
00238  
00239                         reco::BeamSpot tmp_beam = Fit_d_z_likelihood(tmp_par,tmp_error_par);
00240             
00241                         double tmp_par2[7] = {tmp_beam.x0(), tmp_beam.y0(), tmp_beam.z0(),
00242                                                                  tmp_beam.sigmaZ(), tmp_beam.dxdz(), tmp_beam.dydz(),
00243                                                                  tmp_beam.BeamWidthX()};
00244                         
00245                         reco::BeamSpot tmp_lh = Fit_dres_z_likelihood(tmp_par2);
00246 
00247                         if ( isnan(ff_minimum) || std::isinf(ff_minimum) ) {
00248                         
00249                 edm::LogWarning("BSFitter") << "Result is non physical. Log-Likelihood fit did not converge." << std::endl;
00250                                 tmp_lh.setType(reco::BeamSpot::Unknown);
00251                                 return tmp_lh;
00252                         }
00253                         return tmp_lh;
00254                         
00255                 } else {
00256                         
00257                         throw cms::Exception("LogicError")
00258                                 << "Error in BeamSpotProducer/BSFitter: "
00259                                 << "Illegal fit type, options are likelihood or resolution";
00260                 }
00261         } else {
00262 
00263                 throw cms::Exception("LogicError")
00264                         << "Error in BeamSpotProducer/BSFitter: "
00265                         << "Illegal variable type, options are \"z\", \"d\", or \"d*z\"";
00266         }
00267         
00268         
00269 }
00270 
00271 //______________________________________________________________________
00272 reco::BeamSpot BSFitter::Fit_z_likelihood(double *inipar) {
00273 
00274         //std::cout << "Fit_z(double *) called" << std::endl;
00275         //std::cout << "inipar[0]= " << inipar[0] << std::endl;
00276         //std::cout << "inipar[1]= " << inipar[1] << std::endl;
00277         
00278         std::vector<double> par(2,0);
00279         std::vector<double> err(2,0);
00280 
00281         par.push_back(0.0);
00282         par.push_back(7.0);
00283         err.push_back(0.0001);
00284         err.push_back(0.0001);
00285         //par[0] = 0.0; err[0] = 0.0;
00286         //par[1] = 7.0; err[1] = 0.0;
00287 
00288         thePDF->SetPDFs("PDFGauss_z");
00289         thePDF->SetData(fBSvector);
00290         //std::cout << "data loaded"<< std::endl;
00291         
00292         //FunctionMinimum fmin = theFitter->Minimize(*theGausszFcn, par, err, 1, 500, 0.1);
00293         MnUserParameters upar;
00294         upar.Add("X0",    0.,0.);
00295         upar.Add("Y0",    0.,0.);
00296         upar.Add("Z0",    inipar[0],0.001);
00297         upar.Add("sigmaZ",inipar[1],0.001);
00298         
00299         MnMigrad migrad(*thePDF, upar);
00300         
00301         FunctionMinimum fmin = migrad();
00302         ff_minimum = fmin.Fval();
00303         //std::cout << " eval= " << ff_minimum
00304         //                << "/n params[0]= " << fmin.Parameters().Vec()(0) << std::endl;
00305         
00306         /*
00307         TMinuit *gmMinuit = new TMinuit(2); 
00308 
00309         //gmMinuit->SetFCN(z_fcn);
00310         gmMinuit->SetFCN(myFitz_fcn);
00311         
00312         
00313         int ierflg = 0;
00314         double step[2] = {0.001,0.001};
00315         
00316         for (int i = 0; i<2; i++) {   
00317                 gmMinuit->mnparm(i,fpar_name[i].c_str(),inipar[i],step[i],0,0,ierflg);
00318         }
00319         gmMinuit->Migrad();
00320         */
00321         reco::BeamSpot::CovarianceMatrix matrix;
00322 
00323         for (int j = 2 ; j < 4 ; ++j) {
00324                 for(int k = j ; k < 4 ; ++k) {
00325                   matrix(j,k) = fmin.Error().Matrix()(j,k);
00326                 }
00327         }
00328                 
00329         return reco::BeamSpot( reco::BeamSpot::Point(0.,
00330                                                      0.,
00331                                                      fmin.Parameters().Vec()(2)),
00332                                fmin.Parameters().Vec()(3),
00333                                0.,
00334                                0.,
00335                                0.,
00336                                matrix,
00337                                fbeamtype );
00338 }
00339 
00340 //______________________________________________________________________
00341 reco::BeamSpot BSFitter::Fit_z_chi2(double *inipar) {
00342 
00343     // N.B. this fit is not performed anymore but now
00344     // Z is fitted in the same track set used in the d0-phi fit after
00345     // each iteration
00346 
00347     
00348         //std::cout << "Fit_z_chi2() called" << std::endl;
00349         // FIXME: include whole tracker z length for the time being
00350         // ==> add protection and z0 cut
00351         h1z = new TH1F("h1z","z distribution",200,-fMaxZ, fMaxZ);
00352         
00353         std::vector<BSTrkParameters>::const_iterator iparam = fBSvector.begin();
00354 
00355         // HERE check size of track vector
00356         
00357         for( iparam = fBSvector.begin(); iparam != fBSvector.end(); ++iparam) {
00358                 
00359                  h1z->Fill( iparam->z0() );
00360                  //std::cout<<"z0="<<iparam->z0()<<"; sigZ0="<<iparam->sigz0()<<std::endl;
00361         }
00362 
00363         h1z->Fit("gaus","QLM0");
00364         //std::cout << "fitted "<< std::endl;
00365         
00366         TF1 *fgaus = h1z->GetFunction("gaus");
00367         //std::cout << "got function" << std::endl;
00368         double fpar[2] = {fgaus->GetParameter(1), fgaus->GetParameter(2) };
00369         //std::cout<<"Debug fpar[2] = (" <<fpar[0]<<","<<fpar[1]<<")"<<std::endl;
00370         reco::BeamSpot::CovarianceMatrix matrix;
00371         // add matrix values.
00372         matrix(2,2) = fgaus->GetParError(1) * fgaus->GetParError(1);
00373         matrix(3,3) = fgaus->GetParError(2) * fgaus->GetParError(2);
00374         
00375         //delete h1z;
00376 
00377         return reco::BeamSpot( reco::BeamSpot::Point(0.,
00378                                                      0.,
00379                                                      fpar[0]),
00380                                fpar[1],
00381                                0.,
00382                                0.,
00383                                0.,
00384                                matrix,
00385                                fbeamtype );     
00386 
00387         
00388 }
00389 
00390 //______________________________________________________________________
00391 reco::BeamSpot BSFitter::Fit_ited0phi() {
00392 
00393         this->d0phi_Init();
00394     edm::LogInfo("BSFitter") << "number of total input tracks: " << fBSvector.size() << std::endl;
00395         
00396         reco::BeamSpot theanswer;
00397 
00398         if ( (int)fBSvector.size() <= fminNtrks ) {
00399         edm::LogWarning("BSFitter") << "need at least " << fminNtrks << " tracks to run beamline fitter." << std::endl;
00400                 fbeamtype = reco::BeamSpot::Fake;
00401                 theanswer.setType(fbeamtype);
00402                 return theanswer;
00403         }
00404         
00405         theanswer = Fit_d0phi(); //get initial ftmp and ftmprow
00406         if ( goodfit ) fnthite++;
00407         //std::cout << "Initial tempanswer (iteration 0): " << theanswer << std::endl;
00408         
00409         reco::BeamSpot preanswer = theanswer;
00410         
00411         while ( goodfit &&
00412                         ftmprow > fconvergence * fBSvector.size() &&
00413                         ftmprow > fminNtrks  ) {
00414                 
00415                 theanswer = Fit_d0phi();
00416                 fd0cut /= 1.5;
00417                 fchi2cut /= 1.5;
00418                 if ( goodfit &&
00419                         ftmprow > fconvergence * fBSvector.size() &&
00420                         ftmprow > fminNtrks ) {
00421                         preanswer = theanswer;
00422                         //std::cout << "Iteration " << fnthite << ": " << preanswer << std::endl;
00423                         fnthite++;
00424                 }
00425         }
00426         // FIXME: return fit results from previous iteration for both bad fit and for >50% tracks thrown away
00427         //std::cout << "The last iteration, theanswer: " << theanswer << std::endl;
00428         theanswer = preanswer;
00429         //std::cout << "Use previous results from iteration #" << ( fnthite > 0 ? fnthite-1 : 0 ) << std::endl;
00430         //if ( fnthite > 1 ) std::cout << theanswer << std::endl;
00431         
00432     edm::LogInfo("BSFitter") << "Total number of successful iterations = " << ( goodfit ? (fnthite+1) : fnthite ) << std::endl;
00433     if (goodfit) {
00434         fbeamtype = reco::BeamSpot::Tracker;
00435         theanswer.setType(fbeamtype);
00436     }
00437     else {
00438         edm::LogWarning("BSFitter") << "Fit doesn't converge!!!" << std::endl;
00439         fbeamtype = reco::BeamSpot::Unknown;
00440         theanswer.setType(fbeamtype);
00441     }
00442         return theanswer;
00443 }
00444 
00445 
00446 //______________________________________________________________________
00447 reco::BeamSpot BSFitter::Fit_d0phi() {
00448 
00449         //LogDebug ("BSFitter") << " we will use " << fBSvector.size() << " tracks.";
00450     if (fnthite > 0) edm::LogInfo("BSFitter") << " number of tracks used: " << ftmprow << std::endl;
00451         //std::cout << " ftmp = matrix("<<ftmp.GetNrows()<<","<<ftmp.GetNcols()<<")"<<std::endl;
00452         //std::cout << " ftmp(0,0)="<<ftmp(0,0)<<std::endl;
00453         //std::cout << " ftmp(1,0)="<<ftmp(1,0)<<std::endl;
00454         //std::cout << " ftmp(2,0)="<<ftmp(2,0)<<std::endl;
00455         //std::cout << " ftmp(3,0)="<<ftmp(3,0)<<std::endl;
00456         
00457         h1z->Reset();
00458         
00459         
00460         TMatrixD x_result(4,1);
00461         TMatrixDSym V_result(4);
00462         
00463         TMatrixDSym Vint(4);
00464         TMatrixD b(4,1);
00465         
00466         //Double_t weightsum = 0;
00467         
00468         Vint.Zero();
00469         b.Zero();
00470         
00471         TMatrixD g(4,1);
00472         TMatrixDSym temp(4);
00473         
00474         std::vector<BSTrkParameters>::iterator iparam = fBSvector.begin();
00475         ftmprow=0;
00476 
00477         
00478         //edm::LogInfo ("BSFitter") << " test";
00479                 
00480         //std::cout << "BSFitter: fit" << std::endl;
00481         
00482         for( iparam = fBSvector.begin() ;
00483                 iparam != fBSvector.end() ; ++iparam) {
00484                 
00485                 
00486                 //if(i->weight2 == 0) continue;
00487                 
00488                 //if (ftmprow==0) {
00489                 //std::cout << "d0=" << iparam->d0() << " sigd0=" << iparam->sigd0()
00490                 //<< " phi0="<< iparam->phi0() << " z0=" << iparam->z0() << std::endl;
00491                 //std::cout << "d0phi_d0=" << iparam->d0phi_d0() << " d0phi_chi2="<<iparam->d0phi_chi2() << std::endl; 
00492                 //}
00493                 g(0,0) = sin(iparam->phi0());
00494                 g(1,0) = -cos(iparam->phi0());
00495                 g(2,0) = iparam->z0() * g(0,0);
00496                 g(3,0) = iparam->z0() * g(1,0);
00497                 
00498                 
00499                 // average transverse beam width
00500                 double sigmabeam2 = 0.006 * 0.006;
00501                 if (finputBeamWidth > 0 ) sigmabeam2 = finputBeamWidth * finputBeamWidth;
00502         else { edm::LogWarning("BSFitter") << "using in fit beam width = " << sqrt(sigmabeam2) << std::endl; }
00503         
00504                 //double sigma2 = sigmabeam2 +  (iparam->sigd0())* (iparam->sigd0()) / iparam->weight2;
00505                 // this should be 2*sigmabeam2?
00506                 double sigma2 = sigmabeam2 +  (iparam->sigd0())* (iparam->sigd0());
00507 
00508                 TMatrixD ftmptrans(1,4);
00509                 ftmptrans = ftmptrans.Transpose(ftmp);
00510                 TMatrixD dcor = ftmptrans * g;
00511                 double chi2tmp = (iparam->d0() - dcor(0,0)) * (iparam->d0() - dcor(0,0))/sigma2;
00512                 (*iparam) = BSTrkParameters(iparam->z0(),iparam->sigz0(),iparam->d0(),iparam->sigd0(),
00513                                             iparam->phi0(), iparam->pt(),dcor(0,0),chi2tmp);
00514 
00515                 bool pass = true;
00516                 if (fapplyd0cut && fnthite>0 ) {
00517                         if ( std::abs(iparam->d0() - dcor(0,0)) > fd0cut ) pass = false;
00518                         
00519                 }
00520                 if (fapplychi2cut && fnthite>0 ) {
00521                         if ( chi2tmp > fchi2cut ) pass = false;
00522                         
00523                 }
00524                 
00525                 if (pass) {
00526                         temp.Zero();
00527                         for(int j = 0 ; j < 4 ; ++j) {
00528                                 for(int k = j ; k < 4 ; ++k) {
00529                                         temp(j,k) += g(j,0) * g(k,0);
00530                                 }
00531                         }
00532 
00533                 
00534                         Vint += (temp * (1 / sigma2));
00535                         b += (iparam->d0() / sigma2 * g);
00536                         //weightsum += sqrt(i->weight2);
00537                         ftmprow++;
00538             h1z->Fill( iparam->z0() );
00539                 }
00540 
00541                 
00542         }
00543         Double_t determinant;
00544         TDecompBK bk(Vint);
00545         bk.SetTol(1e-11); //FIXME: find a better way to solve x_result
00546         if (!bk.Decompose()) {
00547           goodfit = false;
00548       edm::LogWarning("BSFitter")
00549           << "Decomposition failed, matrix singular ?" << std::endl
00550           << "condition number = " << bk.Condition() << std::endl;
00551         }
00552         else {
00553           V_result = Vint.InvertFast(&determinant);
00554           x_result = V_result  * b;
00555         }
00556         //      for(int j = 0 ; j < 4 ; ++j) {
00557         //        for(int k = 0 ; k < 4 ; ++k) {
00558         //          std::cout<<"V_result("<<j<<","<<k<<")="<<V_result(j,k)<<std::endl;
00559         //        }
00560         //      }
00561         //         for (int j=0;j<4;++j){
00562         //        std::cout<<"x_result("<<j<<",0)="<<x_result(j,0)<<std::endl;
00563         //      }
00564         //LogDebug ("BSFitter") << " d0-phi fit done.";
00565         //std::cout<< " d0-phi fit done." << std::endl;
00566 
00567         h1z->Fit("gaus","QLM0","",h1z->GetMean() -2.*h1z->GetRMS(),h1z->GetMean() +2.*h1z->GetRMS());
00568 
00569         //std::cout << "fitted "<< std::endl;
00570         TF1 *fgaus = h1z->GetFunction("gaus");
00571         //std::cout << "got function" << std::endl;
00572         if (!fgaus){    
00573           edm::LogError("NoBeamSpotFit")<<"gaussian fit failed. no BS d0 fit";          
00574           return reco::BeamSpot();
00575         }
00576         double fpar[2] = {fgaus->GetParameter(1), fgaus->GetParameter(2) };
00577     
00578         reco::BeamSpot::CovarianceMatrix matrix;
00579         // first two parameters
00580         for (int j = 0 ; j < 2 ; ++j) {
00581                 for(int k = j ; k < 2 ; ++k) {
00582                         matrix(j,k) = V_result(j,k);
00583                 }
00584         }
00585         // slope parameters
00586         for (int j = 4 ; j < 6 ; ++j) {
00587                 for(int k = j ; k < 6 ; ++k) {
00588                         matrix(j,k) = V_result(j-2,k-2);
00589                 }
00590         }
00591 
00592     // Z0 and sigmaZ
00593         matrix(2,2) = fgaus->GetParError(1) * fgaus->GetParError(1);
00594         matrix(3,3) = fgaus->GetParError(2) * fgaus->GetParError(2);
00595     
00596         ftmp = x_result;
00597 
00598         // x0 and y0 are *not* x,y at z=0, but actually at z=0
00599         // to correct for this, we need to translate them to z=z0
00600         // using the measured slopes
00601         //
00602         double x0tmp = x_result(0,0);
00603         double y0tmp = x_result(1,0);
00604 
00605         x0tmp += x_result(2,0)*fpar[0];
00606         y0tmp += x_result(3,0)*fpar[0];
00607 
00608 
00609         return reco::BeamSpot( reco::BeamSpot::Point(x0tmp,
00610                                                      y0tmp,
00611                                                      fpar[0]),
00612                            fpar[1],
00613                                x_result(2,0),
00614                                x_result(3,0),
00615                                0.,
00616                                matrix,
00617                                fbeamtype );
00618         
00619 }
00620 
00621 
00622 //______________________________________________________________________
00623 void BSFitter::Setd0Cut_d0phi(double d0cut) {
00624 
00625         fapplyd0cut = true;
00626 
00627         //fBSforCuts = BSfitted;
00628         fd0cut = d0cut;
00629 }
00630 
00631 //______________________________________________________________________
00632 void BSFitter::SetChi2Cut_d0phi(double chi2cut) {
00633 
00634         fapplychi2cut = true;
00635 
00636         //fBSforCuts = BSfitted;
00637         fchi2cut = chi2cut;
00638 }
00639 
00640 //______________________________________________________________________
00641 reco::BeamSpot BSFitter::Fit_d_likelihood(double *inipar) {
00642         
00643 
00644         thePDF->SetPDFs("PDFGauss_d");
00645         thePDF->SetData(fBSvector);
00646 
00647         MnUserParameters upar;
00648         upar.Add("X0",  inipar[0],0.001);
00649         upar.Add("Y0",  inipar[1],0.001);
00650         upar.Add("Z0",    0.,0.001);
00651         upar.Add("sigmaZ",0.,0.001);
00652         upar.Add("dxdz",inipar[2],0.001);
00653         upar.Add("dydz",inipar[3],0.001);
00654                 
00655         
00656         MnMigrad migrad(*thePDF, upar);
00657         
00658         FunctionMinimum fmin = migrad();
00659         ff_minimum = fmin.Fval();
00660 
00661         reco::BeamSpot::CovarianceMatrix matrix;
00662         for (int j = 0 ; j < 6 ; ++j) {
00663                 for(int k = j ; k < 6 ; ++k) {
00664                         matrix(j,k) = fmin.Error().Matrix()(j,k);
00665                 }
00666         }
00667         
00668         return reco::BeamSpot( reco::BeamSpot::Point(fmin.Parameters().Vec()(0),
00669                                                      fmin.Parameters().Vec()(1),
00670                                                      0.),
00671                                0.,
00672                                fmin.Parameters().Vec()(4),
00673                                fmin.Parameters().Vec()(5),
00674                                0.,
00675                                matrix,
00676                                fbeamtype );
00677 }
00678 //______________________________________________________________________
00679 double BSFitter::scanPDF(double *init_pars, int & tracksfixed, int option){
00680 
00681    if(option==1)init_pars[6]=0.0005;  //starting value for any given configuration
00682 
00683    //local vairables with initial values
00684    double fsqrt2pi=0.0;
00685    double d_sig=0.0;
00686    double d_dprime=0.0;
00687    double d_result=0.0;
00688    double z_sig=0.0;
00689    double z_result=0.0;
00690    double function=0.0;
00691    double tot_pdf=0.0;
00692    double last_minvalue=1.0e+10;
00693    double init_bw=-99.99;
00694    int iters=0;
00695 
00696   //used to remove tracks if far away from bs by this
00697    double DeltadCut=0.1000;
00698    if(init_pars[6]<0.0200){DeltadCut=0.0900; } //worked for high 2.36TeV 
00699    if(init_pars[6]<0.0100){DeltadCut=0.0700;}  //just a guesss for 7 TeV but one should scan for actual values
00700 
00701 
00702 std::vector<BSTrkParameters>::const_iterator iparam = fBSvector.begin();
00703 
00704 
00705 if(option==1)iters=500;
00706 if(option==2)iters=1;
00707 
00708 for(int p=0;p<iters;p++){
00709 
00710    if(iters==500)init_pars[6]+=0.0002;
00711     tracksfixed=0;
00712 
00713 for( iparam = fBSvector.begin(); iparam != fBSvector.end(); ++iparam)
00714        {
00715                     fsqrt2pi = sqrt(2.* TMath::Pi());
00716                     d_sig = sqrt(init_pars[6]*init_pars[6] + (iparam->sigd0())*(iparam->sigd0()));
00717                     d_dprime = iparam->d0() - (   (  (init_pars[0] + iparam->z0()*(init_pars[4]))*sin(iparam->phi0()) )
00718                                                - (  (init_pars[1] + iparam->z0()*(init_pars[5]))*cos(iparam->phi0()) ) );
00719 
00720                     //***Remove tracks before the fit which gives low pdf values to blow up the pdf
00721                     if(std::abs(d_dprime)<DeltadCut && option==2){ fBSvectorBW.push_back(*iparam);}
00722 
00723                     d_result = (exp(-(d_dprime*d_dprime)/(2.0*d_sig*d_sig)))/(d_sig*fsqrt2pi);
00724                     z_sig = sqrt(iparam->sigz0() * iparam->sigz0() + init_pars[3]*init_pars[3]);
00725                     z_result = (exp(-((iparam->z0() - init_pars[2])*(iparam->z0() - init_pars[2]))/(2.0*z_sig*z_sig)))/(z_sig*fsqrt2pi);
00726                     tot_pdf=z_result*d_result;
00727 
00728                     //for those trcks which gives problems due to very tiny pdf_d values.
00729                     //Update: This protection will NOT be used with the dprime cut above but still kept here to get
00730                     // the intial value of beam width reasonably
00731                     //A warning will appear if there were any tracks with < 10^-5 for pdf_d so that (d-dprime) cut can be lowered
00732                     if(d_result < 1.0e-05){ tot_pdf=z_result*1.0e-05;
00733                                            //if(option==2)std::cout<<"last Iter  d-d'   =  "<<(std::abs(d_dprime))<<std::endl;
00734                                            tracksfixed++; }
00735 
00736                        function = function + log(tot_pdf);
00737                        tot_pdf=0.0;
00738 
00739 
00740        }//loop over tracks
00741 
00742 
00743        function= -2.0*function;
00744        if(function<last_minvalue){init_bw=init_pars[6];
00745                                   last_minvalue=function; }
00746        function=0.0;
00747    }//loop over beam width
00748 
00749    if(init_bw>0) {
00750     init_bw=init_bw+(0.20*init_bw); //start with 20 % more
00751 
00752    }
00753    else{
00754 
00755        if(option==1){
00756            edm::LogWarning("BSFitter")
00757                <<"scanPDF:====>>>> WARNING***: The initial guess value of Beam width is negative!!!!!!"<<std::endl
00758                <<"scanPDF:====>>>> Assigning beam width a starting value of "<<init_bw<<"  cm"<<std::endl;
00759            init_bw=0.0200;
00760                         
00761        }
00762       }
00763 
00764 
00765     return init_bw;
00766 
00767 }
00768 
00769 //________________________________________________________________________________
00770 reco::BeamSpot BSFitter::Fit_d_z_likelihood(double *inipar, double *error_par) {
00771 
00772       int tracksFailed=0;
00773 
00774       //estimate first guess of beam width and tame 20% extra of it to start
00775       inipar[6]=scanPDF(inipar,tracksFailed,1);
00776       error_par[6]=(inipar[6])*0.20;
00777 
00778 
00779      //Here remove the tracks which give low pdf and fill into a new vector
00780      //std::cout<<"Size of Old vector = "<<(fBSvector.size())<<std::endl;
00781      /* double junk= */ scanPDF(inipar,tracksFailed,2);
00782      //std::cout<<"Size of New vector = "<<(fBSvectorBW.size())<<std::endl;
00783 
00784      //Refill the fBSVector again with new sets of tracks
00785      fBSvector.clear();
00786      std::vector<BSTrkParameters>::const_iterator iparamBW = fBSvectorBW.begin();
00787      for( iparamBW = fBSvectorBW.begin(); iparamBW != fBSvectorBW.end(); ++iparamBW)
00788         {          fBSvector.push_back(*iparamBW); 
00789         }
00790 
00791 
00792         thePDF->SetPDFs("PDFGauss_d*PDFGauss_z");
00793         thePDF->SetData(fBSvector);
00794         MnUserParameters upar;
00795 
00796         upar.Add("X0",  inipar[0],error_par[0]);
00797         upar.Add("Y0",  inipar[1],error_par[1]);
00798         upar.Add("Z0",    inipar[2],error_par[2]);
00799         upar.Add("sigmaZ",inipar[3],error_par[3]);
00800         upar.Add("dxdz",inipar[4],error_par[4]);
00801         upar.Add("dydz",inipar[5],error_par[5]);
00802         upar.Add("BeamWidthX",inipar[6],error_par[6]);
00803 
00804 
00805         MnMigrad migrad(*thePDF, upar);
00806 
00807         FunctionMinimum fmin = migrad();
00808 
00809       // std::cout<<"-----how the fit evoves------"<<std::endl;
00810       // std::cout<<fmin<<std::endl;
00811 
00812         ff_minimum = fmin.Fval();
00813 
00814 
00815         bool ff_nfcn=fmin.HasReachedCallLimit();
00816         bool ff_cov=fmin.HasCovariance();
00817         bool testing=fmin.IsValid();
00818 
00819 
00820         //Print WARNINGS if minimum did not converged
00821         if( ! testing )
00822         {
00823             edm::LogWarning("BSFitter") <<"===========>>>>>** WARNING: MINUIT DID NOT CONVERGES PROPERLY !!!!!!"<<std::endl;
00824             if(ff_nfcn) edm::LogWarning("BSFitter") <<"===========>>>>>** WARNING: No. of Calls Exhausted"<<std::endl;
00825             if(!ff_cov) edm::LogWarning("BSFitter") <<"===========>>>>>** WARNING: Covariance did not found"<<std::endl;
00826         }
00827 
00828         edm::LogInfo("BSFitter") <<"The Total # Tracks used for beam width fit = "<<(fBSvectorBW.size())<<std::endl;
00829 
00830 
00831     //Checks after fit is performed 
00832     double lastIter_pars[7];
00833 
00834    for(int ip=0;ip<7;ip++){ lastIter_pars[ip]=fmin.Parameters().Vec()(ip);
00835                            }
00836 
00837 
00838 
00839     tracksFailed=0;
00840     /* double lastIter_scan= */ scanPDF(lastIter_pars,tracksFailed,2);
00841 
00842    
00843     edm::LogWarning("BSFitter") <<"WARNING: # of tracks which have very low pdf value (pdf_d < 1.0e-05) are  = "<<tracksFailed<<std::endl;
00844 
00845 
00846 
00847         //std::cout << " eval= " << ff_minimum
00848         //                << "/n params[0]= " << fmin.Parameters().Vec()(0) << std::endl;
00849 
00850         reco::BeamSpot::CovarianceMatrix matrix;
00851 
00852         for (int j = 0 ; j < 7 ; ++j) {
00853                 for(int k = j ; k < 7 ; ++k) {
00854                         matrix(j,k) = fmin.Error().Matrix()(j,k);
00855                 }
00856         }
00857 
00858 
00859         return reco::BeamSpot( reco::BeamSpot::Point(fmin.Parameters().Vec()(0),
00860                                                      fmin.Parameters().Vec()(1),
00861                                                      fmin.Parameters().Vec()(2)),
00862                                fmin.Parameters().Vec()(3),
00863                                fmin.Parameters().Vec()(4),
00864                                fmin.Parameters().Vec()(5),
00865                                fmin.Parameters().Vec()(6),
00866                                                               
00867                                matrix,
00868                                fbeamtype );
00869 }
00870 
00871 
00872 //______________________________________________________________________
00873 reco::BeamSpot BSFitter::Fit_dres_z_likelihood(double *inipar) {
00874         
00875         
00876         thePDF->SetPDFs("PDFGauss_d_resolution*PDFGauss_z");
00877         thePDF->SetData(fBSvector);
00878 
00879         MnUserParameters upar;
00880         upar.Add("X0",  inipar[0],0.001);
00881         upar.Add("Y0",  inipar[1],0.001);
00882         upar.Add("Z0",    inipar[2],0.001);
00883         upar.Add("sigmaZ",inipar[3],0.001);
00884         upar.Add("dxdz",inipar[4],0.001);
00885         upar.Add("dydz",inipar[5],0.001);
00886         upar.Add("BeamWidthX",inipar[6],0.0001);
00887         upar.Add("c0",0.0010,0.0001);
00888         upar.Add("c1",0.0090,0.0001);
00889 
00890         // fix beam width
00891         upar.Fix("BeamWidthX");
00892         // number of parameters in fit are 9-1 = 8
00893         
00894         MnMigrad migrad(*thePDF, upar);
00895                 
00896         FunctionMinimum fmin = migrad();
00897         ff_minimum = fmin.Fval();
00898 
00899         reco::BeamSpot::CovarianceMatrix matrix;
00900 
00901         for (int j = 0 ; j < 6 ; ++j) {
00902                 for(int k = j ; k < 6 ; ++k) {
00903                         matrix(j,k) = fmin.Error().Matrix()(j,k);
00904                 }
00905         }
00906 
00907         //std::cout << " fill resolution values" << std::endl;
00908         //std::cout << " matrix size= " << fmin.Error().Matrix().size() << std::endl;
00909         //std::cout << " vec(6)="<< fmin.Parameters().Vec()(6) << std::endl;
00910         //std::cout << " vec(7)="<< fmin.Parameters().Vec()(7) << std::endl;
00911         
00912         fresolution_c0 = fmin.Parameters().Vec()(6);
00913         fresolution_c1 = fmin.Parameters().Vec()(7);
00914         fres_c0_err = sqrt( fmin.Error().Matrix()(6,6) );
00915         fres_c1_err = sqrt( fmin.Error().Matrix()(7,7) );
00916         
00917         for (int j = 6 ; j < 8 ; ++j) {
00918                 for(int k = 6 ; k < 8 ; ++k) {
00919                         fres_matrix(j-6,k-6) = fmin.Error().Matrix()(j,k);
00920                 }
00921         }
00922 
00923         return reco::BeamSpot( reco::BeamSpot::Point(fmin.Parameters().Vec()(0),
00924                                                                          fmin.Parameters().Vec()(1),
00925                                                                          fmin.Parameters().Vec()(2)),
00926                                          fmin.Parameters().Vec()(3),
00927                                          fmin.Parameters().Vec()(4),
00928                                          fmin.Parameters().Vec()(5),
00929                                          inipar[6],
00930                                          matrix,
00931                                          fbeamtype );
00932 }
00933 
00934 
00935
CMSSW_4_4_3_patch1/src/RecoVertex/BeamSpotProducer/src/BSFitter.cc
