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