CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PVFitter.cc
Go to the documentation of this file.
1 
14 #include "RecoVertex/BeamSpotProducer/interface/PVFitter.h"
15 #include "FWCore/ServiceRegistry/interface/Service.h"
16 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
17 #include "CondFormats/BeamSpotObjects/interface/BeamSpotObjects.h"
18 
19 #include "FWCore/Utilities/interface/InputTag.h"
20 #include "DataFormats/Common/interface/View.h"
21 
22 #include "DataFormats/TrackCandidate/interface/TrackCandidate.h"
23 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
24 #include "DataFormats/TrackReco/interface/Track.h"
25 #include "DataFormats/TrackReco/interface/TrackFwd.h"
26 #include "DataFormats/TrackReco/interface/HitPattern.h"
27 #include "DataFormats/VertexReco/interface/Vertex.h"
28 #include "DataFormats/VertexReco/interface/VertexFwd.h"
29 
30 #include "TFitterMinuit.h"
31 #include "Minuit2/FCNBase.h"
32 #include "RecoVertex/BeamSpotProducer/interface/FcnBeamSpotFitPV.h"
33 #include "FWCore/Utilities/interface/isFinite.h"
34 
35 #include "TF1.h"
36 
37 #include "FWCore/MessageLogger/interface/MessageLogger.h"
38 
39 // ----------------------------------------------------------------------
40 // Useful function:
41 // ----------------------------------------------------------------------
42 
43 // static char * formatTime(const std::time_t & t) {
44 // struct std::tm * ptm;
45 // ptm = gmtime(&t);
46 // static char ts[32];
47 // strftime(ts,sizeof(ts),"%Y.%m.%d %H:%M:%S %Z",ptm);
48 // return ts;
49 // }
50 
51 PVFitter::PVFitter(const edm::ParameterSet& iConfig): ftree_(0)
52 {
53 
54  debug_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<bool>("Debug");
55  vertexLabel_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<edm::InputTag>("VertexCollection", edm::InputTag("offlinePrimaryVertices"));
56  do3DFit_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<bool>("Apply3DFit");
57  //writeTxt_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<bool>("WriteAscii");
58  //outputTxt_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<std::string>("AsciiFileName");
59 
60  maxNrVertices_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<unsigned int>("maxNrStoredVertices");
61  minNrVertices_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<unsigned int>("minNrVerticesForFit");
62  minVtxNdf_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("minVertexNdf");
63  maxVtxNormChi2_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("maxVertexNormChi2");
64  minVtxTracks_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<unsigned int>("minVertexNTracks");
65  minVtxWgt_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("minVertexMeanWeight");
66  maxVtxR_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("maxVertexR");
67  maxVtxZ_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("maxVertexZ");
68  errorScale_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("errorScale");
69  sigmaCut_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("nSigmaCut");
70  fFitPerBunchCrossing=iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<bool>("FitPerBunchCrossing");
71 
72  // preset quality cut to "infinite"
73  dynamicQualityCut_ = 1.e30;
74 
75  hPVx = new TH2F("hPVx","PVx vs PVz distribution",200,-maxVtxR_, maxVtxR_, 200, -maxVtxZ_, maxVtxZ_);
76  hPVy = new TH2F("hPVy","PVy vs PVz distribution",200,-maxVtxR_, maxVtxR_, 200, -maxVtxZ_, maxVtxZ_);
77 }
78 
79 PVFitter::~PVFitter() {
80 
81 }
82 
83 
84 void PVFitter::readEvent(const edm::Event& iEvent)
85 {
86 
87 // frun = iEvent.id().run();
88 // const edm::TimeValue_t ftimestamp = iEvent.time().value();
89 // const std::time_t ftmptime = ftimestamp >> 32;
90 
91 // if (fbeginLumiOfFit == -1) freftime[0] = freftime[1] = ftmptime;
92 // if (freftime[0] == 0 || ftmptime < freftime[0]) freftime[0] = ftmptime;
93 // const char* fbeginTime = formatTime(freftime[0]);
94 // sprintf(fbeginTimeOfFit,"%s",fbeginTime);
95 
96 // if (freftime[1] == 0 || ftmptime > freftime[1]) freftime[1] = ftmptime;
97 // const char* fendTime = formatTime(freftime[1]);
98 // sprintf(fendTimeOfFit,"%s",fendTime);
99 
100 // flumi = iEvent.luminosityBlock();
101 // frunFit = frun;
102 
103 // if (fbeginLumiOfFit == -1 || fbeginLumiOfFit > flumi) fbeginLumiOfFit = flumi;
104 // if (fendLumiOfFit == -1 || fendLumiOfFit < flumi) fendLumiOfFit = flumi;
105 // std::cout << "flumi = " <<flumi<<"; fbeginLumiOfFit = " << fbeginLumiOfFit <<"; fendLumiOfFit = "<<fendLumiOfFit<<std::endl;
106 
107  //------ Primary Vertices
108  edm::Handle< reco::VertexCollection > PVCollection;
109  bool hasPVs = false;
110  //edm::View<reco::Vertex> vertices;
111  //const reco::VertexCollection & vertices = 0;
112 
113  if ( iEvent.getByLabel(vertexLabel_, PVCollection ) ) {
114  //pv = *PVCollection;
115  //vertices = *PVCollection;
116  hasPVs = true;
117  }
118  //------
119 
120  if ( hasPVs ) {
121 
122  for (reco::VertexCollection::const_iterator pv = PVCollection->begin(); pv != PVCollection->end(); ++pv ) {
123 
124 
125  //for ( size_t ipv=0; ipv != pv.size(); ++ipv ) {
126 
127  //--- vertex selection
128  if ( pv->isFake() || pv->tracksSize()==0 ) continue;
129  if ( pv->ndof() < minVtxNdf_ || (pv->ndof()+3.)/pv->tracksSize()<2*minVtxWgt_ ) continue;
130  //---
131 
132  hPVx->Fill( pv->x(), pv->z() );
133  hPVy->Fill( pv->y(), pv->z() );
134 
135  //
136  // 3D fit section
137  //
138  // apply additional quality cut
139  if ( pvQuality(*pv)>dynamicQualityCut_ ) continue;
140  // if store exceeds max. size: reduce size and apply new quality cut
141  if ( pvStore_.size()>=maxNrVertices_ ) {
142  compressStore();
143  if ( pvQuality(*pv)>dynamicQualityCut_ ) continue;
144  }
145  //
146  // copy PV to store
147  //
148  int bx = iEvent.bunchCrossing();
149  BeamSpotFitPVData pvData;
150  pvData.bunchCrossing = bx;
151  pvData.position[0] = pv->x();
152  pvData.position[1] = pv->y();
153  pvData.position[2] = pv->z();
154  pvData.posError[0] = pv->xError();
155  pvData.posError[1] = pv->yError();
156  pvData.posError[2] = pv->zError();
157  pvData.posCorr[0] = pv->covariance(0,1)/pv->xError()/pv->yError();
158  pvData.posCorr[1] = pv->covariance(0,2)/pv->xError()/pv->zError();
159  pvData.posCorr[2] = pv->covariance(1,2)/pv->yError()/pv->zError();
160  pvStore_.push_back(pvData);
161 
162  if(ftree_ != 0){
163  theBeamSpotTreeData_.run(iEvent.id().run());
164  theBeamSpotTreeData_.lumi(iEvent.luminosityBlock());
165  theBeamSpotTreeData_.bunchCrossing(bx);
166  theBeamSpotTreeData_.pvData(pvData);
167  ftree_->Fill();
168  }
169 
170  if (fFitPerBunchCrossing) bxMap_[bx].push_back(pvData);
171 
172  }
173 
174  }
175 
176 
177 
178 
179 }
180 
181 void PVFitter::setTree(TTree* tree){
182  ftree_ = tree;
183  theBeamSpotTreeData_.branch(ftree_);
184 }
185 
186 bool PVFitter::runBXFitter() {
187 
188  edm::LogInfo("PVFitter") << " Number of bunch crossings: " << bxMap_.size() << std::endl;
189 
190  bool fit_ok = true;
191 
192  for ( std::map<int,std::vector<BeamSpotFitPVData> >::const_iterator pvStore = bxMap_.begin();
193  pvStore!=bxMap_.end(); ++pvStore) {
194 
195  // first set null beam spot in case
196  // fit fails
197  fbspotMap[pvStore->first] = reco::BeamSpot();
198 
199  edm::LogInfo("PVFitter") << " Number of PVs collected for PVFitter: " << (pvStore->second).size() << " in bx: " << pvStore->first << std::endl;
200 
201  if ( (pvStore->second).size() <= minNrVertices_ ) {
202  edm::LogWarning("PVFitter") << " not enough PVs, continue" << std::endl;
203  fit_ok = false;
204  continue;
205  }
206 
207  //bool fit_ok = false;
208  edm::LogInfo("PVFitter") << "Calculating beam spot with PVs ..." << std::endl;
209 
210  //
211  // LL function and fitter
212  //
213  FcnBeamSpotFitPV* fcn = new FcnBeamSpotFitPV(pvStore->second);
214  TFitterMinuit minuitx;
215  minuitx.SetMinuitFCN(fcn);
216  //
217  // fit parameters: positions, widths, x-y correlations, tilts in xz and yz
218  //
219  minuitx.SetParameter(0,"x",0.,0.02,-10.,10.);
220  minuitx.SetParameter(1,"y",0.,0.02,-10.,10.);
221  minuitx.SetParameter(2,"z",0.,0.20,-30.,30.);
222  minuitx.SetParameter(3,"ex",0.015,0.01,0.,10.);
223  minuitx.SetParameter(4,"corrxy",0.,0.02,-1.,1.);
224  minuitx.SetParameter(5,"ey",0.015,0.01,0.,10.);
225  minuitx.SetParameter(6,"dxdz",0.,0.0002,-0.1,0.1);
226  minuitx.SetParameter(7,"dydz",0.,0.0002,-0.1,0.1);
227  minuitx.SetParameter(8,"ez",1.,0.1,0.,30.);
228  minuitx.SetParameter(9,"scale",errorScale_,errorScale_/10.,errorScale_/2.,errorScale_*2.);
229  //
230  // first iteration without correlations
231  //
232  int ierr(0);
233  minuitx.FixParameter(4);
234  minuitx.FixParameter(6);
235  minuitx.FixParameter(7);
236  minuitx.FixParameter(9);
237  minuitx.SetMaxIterations(100);
238  // minuitx.SetPrintLevel(3);
239  minuitx.SetPrintLevel(0);
240  minuitx.CreateMinimizer();
241  ierr = minuitx.Minimize();
242  if ( ierr ) {
243  edm::LogInfo("PVFitter") << "3D beam spot fit failed in 1st iteration" << std::endl;
244  fit_ok = false;
245  continue;
246  }
247  //
248  // refit with harder selection on vertices
249  //
250  fcn->setLimits(minuitx.GetParameter(0)-sigmaCut_*minuitx.GetParameter(3),
251  minuitx.GetParameter(0)+sigmaCut_*minuitx.GetParameter(3),
252  minuitx.GetParameter(1)-sigmaCut_*minuitx.GetParameter(5),
253  minuitx.GetParameter(1)+sigmaCut_*minuitx.GetParameter(5),
254  minuitx.GetParameter(2)-sigmaCut_*minuitx.GetParameter(8),
255  minuitx.GetParameter(2)+sigmaCut_*minuitx.GetParameter(8));
256  ierr = minuitx.Minimize();
257  if ( ierr ) {
258  edm::LogInfo("PVFitter") << "3D beam spot fit failed in 2nd iteration" << std::endl;
259  fit_ok = false;
260  continue;
261  }
262  //
263  // refit with correlations
264  //
265  minuitx.ReleaseParameter(4);
266  minuitx.ReleaseParameter(6);
267  minuitx.ReleaseParameter(7);
268 
269  ierr = minuitx.Minimize();
270  if ( ierr ) {
271  edm::LogInfo("PVFitter") << "3D beam spot fit failed in 3rd iteration" << std::endl;
272  fit_ok = false;
273  continue;
274  }
275  // refit with floating scale factor
276  // minuitx.ReleaseParameter(9);
277  // minuitx.Minimize();
278 
279  //minuitx.PrintResults(0,0);
280 
281  fwidthX = minuitx.GetParameter(3);
282  fwidthY = minuitx.GetParameter(5);
283  fwidthZ = minuitx.GetParameter(8);
284  fwidthXerr = minuitx.GetParError(3);
285  fwidthYerr = minuitx.GetParError(5);
286  fwidthZerr = minuitx.GetParError(8);
287 
288  reco::BeamSpot::CovarianceMatrix matrix;
289  // need to get the full cov matrix
290  matrix(0,0) = pow( minuitx.GetParError(0), 2);
291  matrix(1,1) = pow( minuitx.GetParError(1), 2);
292  matrix(2,2) = pow( minuitx.GetParError(2), 2);
293  matrix(3,3) = fwidthZerr * fwidthZerr;
294  matrix(4,4) = pow( minuitx.GetParError(6), 2);
295  matrix(5,5) = pow( minuitx.GetParError(7), 2);
296  matrix(6,6) = fwidthXerr * fwidthXerr;
297 
298  fbeamspot = reco::BeamSpot( reco::BeamSpot::Point(minuitx.GetParameter(0),
299  minuitx.GetParameter(1),
300  minuitx.GetParameter(2) ),
301  fwidthZ,
302  minuitx.GetParameter(6), minuitx.GetParameter(7),
303  fwidthX,
304  matrix );
305  fbeamspot.setBeamWidthX( fwidthX );
306  fbeamspot.setBeamWidthY( fwidthY );
307  fbeamspot.setType( reco::BeamSpot::Tracker );
308 
309  fbspotMap[pvStore->first] = fbeamspot;
310  edm::LogInfo("PVFitter") << "3D PV fit done for this bunch crossing."<<std::endl;
311  minuitx.Clear();
312  //delete fcn;
313  fit_ok = fit_ok & true;
314  }
315 
316  return fit_ok;
317 }
318 
319 
320 bool PVFitter::runFitter() {
321 
322  edm::LogInfo("PVFitter") << " Number of PVs collected for PVFitter: " << pvStore_.size() << std::endl;
323 
324  if ( pvStore_.size() <= minNrVertices_ ) return false;
325 
326  //bool fit_ok = false;
327 
328  if ( ! do3DFit_ ) {
329  TH1F *h1PVx = (TH1F*) hPVx->ProjectionX("h1PVx", 0, -1, "e");
330  TH1F *h1PVy = (TH1F*) hPVy->ProjectionX("h1PVy", 0, -1, "e");
331  TH1F *h1PVz = (TH1F*) hPVx->ProjectionY("h1PVz", 0, -1, "e");
332 
333  h1PVx->Fit("gaus","QLM0");
334  h1PVy->Fit("gaus","QLM0");
335  h1PVz->Fit("gaus","QLM0");
336 
337  TF1 *gausx = h1PVx->GetFunction("gaus");
338  TF1 *gausy = h1PVy->GetFunction("gaus");
339  TF1 *gausz = h1PVz->GetFunction("gaus");
340 
341  fwidthX = gausx->GetParameter(2);
342  fwidthY = gausy->GetParameter(2);
343  fwidthZ = gausz->GetParameter(2);
344  fwidthXerr = gausx->GetParError(2);
345  fwidthYerr = gausy->GetParError(2);
346  fwidthZerr = gausz->GetParError(2);
347 
348  reco::BeamSpot::CovarianceMatrix matrix;
349  matrix(2,2) = gausz->GetParError(1) * gausz->GetParError(1);
350  matrix(3,3) = fwidthZerr * fwidthZerr;
351  matrix(6,6) = fwidthXerr * fwidthXerr;
352 
353  fbeamspot = reco::BeamSpot( reco::BeamSpot::Point(gausx->GetParameter(1),
354  gausy->GetParameter(1),
355  gausz->GetParameter(1) ),
356  fwidthZ,
357  0., 0.,
358  fwidthX,
359  matrix );
360  fbeamspot.setBeamWidthX( fwidthX );
361  fbeamspot.setBeamWidthY( fwidthY );
362  fbeamspot.setType(reco::BeamSpot::Tracker);
363 
364  }
365  else { // do 3D fit
366  //
367  // LL function and fitter
368  //
369  FcnBeamSpotFitPV* fcn = new FcnBeamSpotFitPV(pvStore_);
370  TFitterMinuit minuitx;
371  minuitx.SetMinuitFCN(fcn);
372  //
373  // fit parameters: positions, widths, x-y correlations, tilts in xz and yz
374  //
375  minuitx.SetParameter(0,"x",0.,0.02,-10.,10.);
376  minuitx.SetParameter(1,"y",0.,0.02,-10.,10.);
377  minuitx.SetParameter(2,"z",0.,0.20,-30.,30.);
378  //minuitx.SetParameter(3,"ex",0.015,0.01,0.,10.);
379  //minuitx.SetParameter(3,"ex",0.015,0.01,0.0001,10.);
380  minuitx.SetParameter(3,"ex",0.005,0.0005,0.0001,0.05);
381  minuitx.SetParameter(4,"corrxy",0.,0.02,-1.,1.);
382  //minuitx.SetParameter(5,"ey",0.015,0.01,0.,10.);
383  //minuitx.SetParameter(5,"ey",0.015,0.01,0.0001,10.);
384  minuitx.SetParameter(5,"ey",0.005,0.0005,0.0001,0.05);
385  minuitx.SetParameter(6,"dxdz",0.,0.0002,-0.1,0.1);
386  minuitx.SetParameter(7,"dydz",0.,0.0002,-0.1,0.1);
387  //minuitx.SetParameter(8,"ez",1.,0.1,0.,30.);
388  minuitx.SetParameter(8,"ez",1.,0.1,1.0,30.);
389  minuitx.SetParameter(9,"scale",errorScale_,errorScale_/10.,errorScale_/2.,errorScale_*2.);
390  //
391  // first iteration without correlations
392  //
393  int ierr(0);
394  minuitx.FixParameter(4);
395  minuitx.FixParameter(6);
396  minuitx.FixParameter(7);
397  minuitx.FixParameter(9);
398  minuitx.SetMaxIterations(100);
399 // minuitx.SetPrintLevel(3);
400  minuitx.SetPrintLevel(0);
401  minuitx.CreateMinimizer();
402  ierr = minuitx.Minimize();
403  if ( ierr ) {
404  edm::LogWarning("PVFitter") << "3D beam spot fit failed in 1st iteration" << std::endl;
405  return false;
406  }
407  //
408  // refit with harder selection on vertices
409  //
410  fcn->setLimits(minuitx.GetParameter(0)-sigmaCut_*minuitx.GetParameter(3),
411  minuitx.GetParameter(0)+sigmaCut_*minuitx.GetParameter(3),
412  minuitx.GetParameter(1)-sigmaCut_*minuitx.GetParameter(5),
413  minuitx.GetParameter(1)+sigmaCut_*minuitx.GetParameter(5),
414  minuitx.GetParameter(2)-sigmaCut_*minuitx.GetParameter(8),
415  minuitx.GetParameter(2)+sigmaCut_*minuitx.GetParameter(8));
416  ierr = minuitx.Minimize();
417  if ( ierr ) {
418  edm::LogWarning("PVFitter") << "3D beam spot fit failed in 2nd iteration" << std::endl;
419  return false;
420  }
421  //
422  // refit with correlations
423  //
424  minuitx.ReleaseParameter(4);
425  minuitx.ReleaseParameter(6);
426  minuitx.ReleaseParameter(7);
427  ierr = minuitx.Minimize();
428  if ( ierr ) {
429  edm::LogWarning("PVFitter") << "3D beam spot fit failed in 3rd iteration" << std::endl;
430  return false;
431  }
432  // refit with floating scale factor
433  // minuitx.ReleaseParameter(9);
434  // minuitx.Minimize();
435 
436  //minuitx.PrintResults(0,0);
437 
438  fwidthX = minuitx.GetParameter(3);
439  fwidthY = minuitx.GetParameter(5);
440  fwidthZ = minuitx.GetParameter(8);
441  fwidthXerr = minuitx.GetParError(3);
442  fwidthYerr = minuitx.GetParError(5);
443  fwidthZerr = minuitx.GetParError(8);
444 
445  // check errors on widths and sigmaZ for nan
446  if ( edm::isNotFinite(fwidthXerr) || edm::isNotFinite(fwidthYerr) || edm::isNotFinite(fwidthZerr) ) {
447  edm::LogWarning("PVFitter") << "3D beam spot fit returns nan in 3rd iteration" << std::endl;
448  return false;
449  }
450 
451  reco::BeamSpot::CovarianceMatrix matrix;
452  // need to get the full cov matrix
453  matrix(0,0) = pow( minuitx.GetParError(0), 2);
454  matrix(1,1) = pow( minuitx.GetParError(1), 2);
455  matrix(2,2) = pow( minuitx.GetParError(2), 2);
456  matrix(3,3) = fwidthZerr * fwidthZerr;
457  matrix(6,6) = fwidthXerr * fwidthXerr;
458 
459  fbeamspot = reco::BeamSpot( reco::BeamSpot::Point(minuitx.GetParameter(0),
460  minuitx.GetParameter(1),
461  minuitx.GetParameter(2) ),
462  fwidthZ,
463  minuitx.GetParameter(6), minuitx.GetParameter(7),
464  fwidthX,
465  matrix );
466  fbeamspot.setBeamWidthX( fwidthX );
467  fbeamspot.setBeamWidthY( fwidthY );
468  fbeamspot.setType(reco::BeamSpot::Tracker);
469  }
470 
471  return true; //FIXME: Need to add quality test for the fit results!
472 }
473 
474 void PVFitter::dumpTxtFile(){
475 /*
476  fasciiFile << "Runnumber " << frun << std::endl;
477  fasciiFile << "BeginTimeOfFit " << fbeginTimeOfFit << std::endl;
478  fasciiFile << "EndTimeOfFit " << fendTimeOfFit << std::endl;
479  fasciiFile << "LumiRange " << fbeginLumiOfFit << " - " << fendLumiOfFit << std::endl;
480  fasciiFile << "Type " << fbeamspot.type() << std::endl;
481  fasciiFile << "X0 " << fbeamspot.x0() << std::endl;
482  fasciiFile << "Y0 " << fbeamspot.y0() << std::endl;
483  fasciiFile << "Z0 " << fbeamspot.z0() << std::endl;
484  fasciiFile << "sigmaZ0 " << fbeamspot.sigmaZ() << std::endl;
485  fasciiFile << "dxdz " << fbeamspot.dxdz() << std::endl;
486  fasciiFile << "dydz " << fbeamspot.dydz() << std::endl;
487  if (inputBeamWidth_ > 0 ) {
488  fasciiFile << "BeamWidthX " << inputBeamWidth_ << std::endl;
489  fasciiFile << "BeamWidthY " << inputBeamWidth_ << std::endl;
490  } else {
491  fasciiFile << "BeamWidthX " << fbeamspot.BeamWidthX() << std::endl;
492  fasciiFile << "BeamWidthY " << fbeamspot.BeamWidthY() << std::endl;
493  }
494 
495  for (int i = 0; i<6; ++i) {
496  fasciiFile << "Cov("<<i<<",j) ";
497  for (int j=0; j<7; ++j) {
498  fasciiFile << fbeamspot.covariance(i,j) << " ";
499  }
500  fasciiFile << std::endl;
501  }
502  // beam width error
503  if (inputBeamWidth_ > 0 ) {
504  fasciiFile << "Cov(6,j) 0 0 0 0 0 0 " << "1e-4" << std::endl;
505  } else {
506  fasciiFile << "Cov(6,j) 0 0 0 0 0 0 " << fbeamspot.covariance(6,6) << std::endl;
507  }
508  fasciiFile << "EmittanceX " << fbeamspot.emittanceX() << std::endl;
509  fasciiFile << "EmittanceY " << fbeamspot.emittanceY() << std::endl;
510  fasciiFile << "BetaStar " << fbeamspot.betaStar() << std::endl;
511 
512 */
513 }
514 
515 
516 void
517 PVFitter::compressStore ()
518 {
519  //
520  // fill vertex qualities
521  //
522  pvQualities_.resize(pvStore_.size());
523  for ( unsigned int i=0; i<pvStore_.size(); ++i ) pvQualities_[i] = pvQuality(pvStore_[i]);
524  sort(pvQualities_.begin(),pvQualities_.end());
525  //
526  // Set new quality cut to median. This cut will be used to reduce the
527  // number of vertices in the store and also apply to all new vertices
528  // until the next reset
529  //
530  dynamicQualityCut_ = pvQualities_[pvQualities_.size()/2];
531  //
532  // remove all vertices failing the cut from the store
533  // (to be moved to a more efficient memory management!)
534  //
535  unsigned int iwrite(0);
536  for ( unsigned int i=0; i<pvStore_.size(); ++i ) {
537  if ( pvQuality(pvStore_[i])>dynamicQualityCut_ ) continue;
538  if ( i!=iwrite ) pvStore_[iwrite] = pvStore_[i];
539  ++iwrite;
540  }
541  pvStore_.resize(iwrite);
542  edm::LogInfo("PVFitter") << "Reduced primary vertex store size to "
543  << pvStore_.size() << " ; new dynamic quality cut = "
544  << dynamicQualityCut_ << std::endl;
545 
546 }
547 
548 double
549 PVFitter::pvQuality (const reco::Vertex& pv) const
550 {
551  //
552  // determinant of the transverse part of the PV covariance matrix
553  //
554  return
555  pv.covariance(0,0)*pv.covariance(1,1)-
556  pv.covariance(0,1)*pv.covariance(0,1);
557 }
558 
559 double
560 PVFitter::pvQuality (const BeamSpotFitPVData& pv) const
561 {
562  //
563  // determinant of the transverse part of the PV covariance matrix
564  //
565  double ex = pv.posError[0];
566  double ey = pv.posError[1];
567  return ex*ex*ey*ey*(1-pv.posCorr[0]*pv.posCorr[0]);
568 }
569 
edm::EventID::run
RunNumber_t run() const
Definition: EventID.h:42
PVFitter::pvQualities_
std::vector< double > pvQualities_
Definition: PVFitter.h:199
reco::BeamSpot::CovarianceMatrix
math::Error< dimension >::type CovarianceMatrix
Definition: BeamSpot.h:32
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
i
int i
Definition: DBlmapReader.cc:9
PVFitter::debug_
bool debug_
Definition: PVFitter.h:139
PVFitter::maxVtxNormChi2_
double maxVtxNormChi2_
Definition: PVFitter.h:148
PVFitter::setTree
void setTree(TTree *tree)
Definition: PVFitter.cc:181
PVFitter::bxMap_
std::map< int, std::vector< BeamSpotFitPVData > > bxMap_
Definition: PVFitter.h:197
BeamSpotTreeData::run
void run(unsigned int run)
Definition: BeamSpotTreeData.h:17
PVFitter::vertexLabel_
edm::InputTag vertexLabel_
Definition: PVFitter.h:141
PVFitter::compressStore
void compressStore()
reduce size of primary vertex cache by increasing quality limit
Definition: PVFitter.cc:517
PVFitter::minVtxWgt_
double minVtxWgt_
Definition: PVFitter.h:150
PVFitter::minVtxNdf_
double minVtxNdf_
Definition: PVFitter.h:147
edm::EventBase::bunchCrossing
int bunchCrossing() const
Definition: EventBase.h:62
edm::EventBase::luminosityBlock
edm::LuminosityBlockNumber_t luminosityBlock() const
Definition: EventBase.h:59
reco::BeamSpot::Point
math::XYZPoint Point
point in the space
Definition: BeamSpot.h:30
reco::Vertex::covariance
double covariance(int i, int j) const
(i, j)-th element of error matrix, i, j = 0, ... 2
Definition: Vertex.h:110
PVFitter::fwidthY
double fwidthY
Definition: PVFitter.h:177
PVFitter::errorScale_
double errorScale_
Definition: PVFitter.h:153
PVFitter::do3DFit_
bool do3DFit_
Definition: PVFitter.h:140
PVFitter::dynamicQualityCut_
double dynamicQualityCut_
Definition: PVFitter.h:198
PVFitter::fwidthYerr
double fwidthYerr
Definition: PVFitter.h:180
PVFitter::fwidthZ
double fwidthZ
Definition: PVFitter.h:178
PVFitter::fwidthXerr
double fwidthXerr
Definition: PVFitter.h:179
Association.map
dictionary map
Definition: Association.py:205
reco::BeamSpot::setType
void setType(BeamType type)
set beam type
Definition: BeamSpot.h:132
BeamSpotTreeData::bunchCrossing
void bunchCrossing(unsigned int bunchCrossing)
Definition: BeamSpotTreeData.h:19
PVFitter::fbeamspot
reco::BeamSpot fbeamspot
Definition: PVFitter.h:133
BeamSpotTreeData::lumi
void lumi(unsigned int lumi)
Definition: BeamSpotTreeData.h:18
iEvent
int iEvent
Definition: GenABIO.cc:243
PVFitter::pvStore_
std::vector< BeamSpotFitPVData > pvStore_
Definition: PVFitter.h:196
PVFitter::readEvent
void readEvent(const edm::Event &iEvent)
Definition: PVFitter.cc:84
edm::isNotFinite
bool isNotFinite(T x)
Definition: isFinite.h:10
reco::BeamSpot::setBeamWidthY
void setBeamWidthY(double v)
Definition: BeamSpot.h:110
PVFitter::theBeamSpotTreeData_
BeamSpotTreeData theBeamSpotTreeData_
Definition: PVFitter.h:201
PVFitter::fFitPerBunchCrossing
bool fFitPerBunchCrossing
Definition: PVFitter.h:135
PVFitter::hPVy
TH2F * hPVy
Definition: PVFitter.h:160
BeamSpotTreeData::pvData
void pvData(const BeamSpotFitPVData &pvData)
Definition: BeamSpotTreeData.h:20
PVFitter::ftree_
TTree * ftree_
Definition: PVFitter.h:162
PVFitter::fbspotMap
std::map< int, reco::BeamSpot > fbspotMap
Definition: PVFitter.h:134
PVFitter::pvQuality
double pvQuality(const reco::Vertex &pv) const
vertex quality measure
Definition: PVFitter.cc:549
PVFitter::maxVtxR_
double maxVtxR_
Definition: PVFitter.h:151
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:361
PVFitter::maxVtxZ_
double maxVtxZ_
Definition: PVFitter.h:152
PVFitter::maxNrVertices_
unsigned int maxNrVertices_
Definition: PVFitter.h:145
PVFitter::fwidthZerr
double fwidthZerr
Definition: PVFitter.h:181
BeamSpotTreeData::branch
void branch(TTree *tree)
Definition: BeamSpotTreeData.cc:10
PVFitter::PVFitter
PVFitter()
Definition: PVFitter.h:41
PVFitter::minVtxTracks_
unsigned int minVtxTracks_
Definition: PVFitter.h:149
fcn
void fcn(int &, double *, double &, double *, int)
Definition: LASBarrelAlgorithm.cc:343
PVFitter::runFitter
bool runFitter()
Definition: PVFitter.cc:320
PVFitter::hPVx
TH2F * hPVx
Definition: PVFitter.h:160
PVFitter::minNrVertices_
unsigned int minNrVertices_
Definition: PVFitter.h:146
FcnBeamSpotFitPV::setLimits
void setLimits(float xmin, float xmax, float ymin, float ymax, float zmin, float zmax)
Definition: FcnBeamSpotFitPV.cc:22
PVFitter::runBXFitter
bool runBXFitter()
Definition: PVFitter.cc:186
PVFitter::sigmaCut_
double sigmaCut_
Definition: PVFitter.h:154
edm::EventBase::id
edm::EventID id() const
Definition: EventBase.h:56
reco::BeamSpot::setBeamWidthX
void setBeamWidthX(double v)
Definition: BeamSpot.h:109
PVFitter::~PVFitter
virtual ~PVFitter()
Definition: PVFitter.cc:79
PVFitter::fwidthX
double fwidthX
Definition: PVFitter.h:176
PVFitter::dumpTxtFile
void dumpTxtFile()
Definition: PVFitter.cc:474
findQualityFiles.size
tuple size
Write out results.
Definition: findQualityFiles.py:442
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
Generated for CMSSW Reference Manual by doxygen 1.8.5