PVFitter Class Reference

#include <PVFitter.h>

Public Member Functions
void	compressStore ()
	reduce size of primary vertex cache by increasing quality limit More...
void	dumpTxtFile ()
void	FitPerBunchCrossing ()
reco::BeamSpot	getBeamSpot ()
std::map< int, reco::BeamSpot >	getBeamSpotMap ()
int *	getFitLSRange ()
int	getNPVs ()
const std::map< int, int > &	getNPVsperBX ()
std::vector< BeamSpotFitPVData >	getpvStore ()
time_t *	getRefTime ()
double	getWidthX ()
double	getWidthXerr ()
double	getWidthY ()
double	getWidthYerr ()
double	getWidthZ ()
double	getWidthZerr ()
bool	IsFitPerBunchCrossing ()
	PVFitter ()
	PVFitter (const edm::ParameterSet &iConfig)
double	pvQuality (const reco::Vertex &pv) const
	vertex quality measure More...
double	pvQuality (const BeamSpotFitPVData &pv) const
	vertex quality measure More...
void	readEvent (const edm::Event &iEvent)
void	resetAll ()
void	resetLSRange ()
void	resetRefTime ()
void	resizepvStore (unsigned int rmSize)
bool	runBXFitter ()
bool	runFitter ()
void	setFitLSRange (int ls0, int ls1)
void	setRefTime (std::time_t t0, std::time_t t1)
void	setTree (TTree *tree)
virtual	~PVFitter ()

Private Attributes
std::map< int, std::vector < BeamSpotFitPVData > >	bxMap_
bool	debug_
bool	do3DFit_
double	dynamicQualityCut_
double	errorScale_
std::ofstream	fasciiFile
reco::BeamSpot	fbeamspot
int	fbeginLumiOfFit
char	fbeginTimeOfFit [32]
std::map< int, reco::BeamSpot >	fbspotMap
double	fdxdz
double	fdxdzErr
double	fdydz
double	fdydzErr
int	fendLumiOfFit
char	fendTimeOfFit [32]
bool	fFitPerBunchCrossing
int	flumi
std::time_t	freftime [2]
int	frun
int	frunFit
double	fsigmaZ
double	fsigmaZErr
TTree *	ftree_
double	fwidthX
double	fwidthXerr
double	fwidthY
double	fwidthYerr
double	fwidthZ
double	fwidthZerr
double	fx
double	fxErr
double	fy
double	fyErr
double	fz
double	fzErr
TH2F *	hPVx
TH2F *	hPVy
unsigned int	maxNrVertices_
double	maxVtxNormChi2_
double	maxVtxR_
double	maxVtxZ_
unsigned int	minNrVertices_
double	minVtxNdf_
unsigned int	minVtxTracks_
double	minVtxWgt_
std::map< int, int >	npvsmap_
std::string	outputTxt_
std::vector< double >	pvQualities_
std::vector< BeamSpotFitPVData >	pvStore_
double	sigmaCut_
BeamSpotTreeData	theBeamSpotTreeData_
edm::InputTag	vertexLabel_
bool	writeTxt_

Detailed Description

Definition at line 39 of file PVFitter.h.

Constructor & Destructor Documentation

PVFitter::PVFitter ( )

inline

Definition at line 41 of file PVFitter.h.

{}

PVFitter::PVFitter

(

const edm::ParameterSet &

iConfig

)

---

class: PVFitter.cc package: RecoVertex/BeamSpotProducer

author: Francisco Yumiceva, Fermilab (yumic.nosp@m.eva@.nosp@m.fnal..nosp@m.gov) Geng-Yuan Jeng, UC Riverside (Geng-.nosp@m.Yuan.nosp@m..Jeng.nosp@m.@cer.nosp@m.n.ch)

version

Id:

PVFitter.cc,v 1.19 2011/12/05 23:14:31 burkett Exp

---

Definition at line 50 of file PVFitter.cc.

References debug_, do3DFit_, dynamicQualityCut_, errorScale_, fFitPerBunchCrossing, edm::ParameterSet::getParameter(), hPVx, hPVy, maxNrVertices_, maxVtxNormChi2_, maxVtxR_, maxVtxZ_, minNrVertices_, minVtxNdf_, minVtxTracks_, minVtxWgt_, sigmaCut_, and vertexLabel_.

                                                : ftree_(0)
{

  debug_             = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<bool>("Debug");
  vertexLabel_     = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<edm::InputTag>("VertexCollection", edm::InputTag("offlinePrimaryVertices"));
  do3DFit_           = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<bool>("Apply3DFit");
  //writeTxt_          = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<bool>("WriteAscii");
  //outputTxt_         = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<std::string>("AsciiFileName");

  maxNrVertices_     = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<unsigned int>("maxNrStoredVertices");
  minNrVertices_     = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<unsigned int>("minNrVerticesForFit");
  minVtxNdf_         = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("minVertexNdf");
  maxVtxNormChi2_    = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("maxVertexNormChi2");
  minVtxTracks_      = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<unsigned int>("minVertexNTracks");
  minVtxWgt_         = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("minVertexMeanWeight");
  maxVtxR_           = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("maxVertexR");
  maxVtxZ_           = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("maxVertexZ");
  errorScale_        = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("errorScale");
  sigmaCut_          = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("nSigmaCut");
  fFitPerBunchCrossing=iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<bool>("FitPerBunchCrossing");

  // preset quality cut to "infinite"
  dynamicQualityCut_ = 1.e30;

  hPVx = new TH2F("hPVx","PVx vs PVz distribution",200,-maxVtxR_, maxVtxR_, 200, -maxVtxZ_, maxVtxZ_);
  hPVy = new TH2F("hPVy","PVy vs PVz distribution",200,-maxVtxR_, maxVtxR_, 200, -maxVtxZ_, maxVtxZ_);
}

PVFitter::~PVFitter ( )

virtual

Definition at line 78 of file PVFitter.cc.

                    {

}

Member Function Documentation

void PVFitter::compressStore

(

)

reduce size of primary vertex cache by increasing quality limit

Definition at line 514 of file PVFitter.cc.

References dynamicQualityCut_, i, pvQualities_, pvQuality(), pvStore_, and python.multivaluedict::sort().

Referenced by readEvent().

{
  //
  // fill vertex qualities
  //
  pvQualities_.resize(pvStore_.size());
  for ( unsigned int i=0; i<pvStore_.size(); ++i )  pvQualities_[i] = pvQuality(pvStore_[i]);
  sort(pvQualities_.begin(),pvQualities_.end());
  //
  // Set new quality cut to median. This cut will be used to reduce the
  // number of vertices in the store and also apply to all new vertices
  // until the next reset
  //
  dynamicQualityCut_ = pvQualities_[pvQualities_.size()/2];
  //
  // remove all vertices failing the cut from the store
  //   (to be moved to a more efficient memory management!)
  //
  unsigned int iwrite(0);
  for ( unsigned int i=0; i<pvStore_.size(); ++i ) {
    if ( pvQuality(pvStore_[i])>dynamicQualityCut_ )  continue;
    if ( i!=iwrite )  pvStore_[iwrite] = pvStore_[i];
    ++iwrite;
  }
  pvStore_.resize(iwrite);
  edm::LogInfo("PVFitter") << "Reduced primary vertex store size to "
                           << pvStore_.size() << " ; new dynamic quality cut = "
                           << dynamicQualityCut_ << std::endl;

}

void PVFitter::dumpTxtFile

(

)

Definition at line 471 of file PVFitter.cc.

                          {
/*
  fasciiFile << "Runnumber " << frun << std::endl;
  fasciiFile << "BeginTimeOfFit " << fbeginTimeOfFit << std::endl;
  fasciiFile << "EndTimeOfFit " << fendTimeOfFit << std::endl;
  fasciiFile << "LumiRange " << fbeginLumiOfFit << " - " << fendLumiOfFit << std::endl;
  fasciiFile << "Type " << fbeamspot.type() << std::endl;
  fasciiFile << "X0 " << fbeamspot.x0() << std::endl;
  fasciiFile << "Y0 " << fbeamspot.y0() << std::endl;
  fasciiFile << "Z0 " << fbeamspot.z0() << std::endl;
  fasciiFile << "sigmaZ0 " << fbeamspot.sigmaZ() << std::endl;
  fasciiFile << "dxdz " << fbeamspot.dxdz() << std::endl;
  fasciiFile << "dydz " << fbeamspot.dydz() << std::endl;
  if (inputBeamWidth_ > 0 ) {
    fasciiFile << "BeamWidthX " << inputBeamWidth_ << std::endl;
    fasciiFile << "BeamWidthY " << inputBeamWidth_ << std::endl;
  } else {
    fasciiFile << "BeamWidthX " << fbeamspot.BeamWidthX() << std::endl;
    fasciiFile << "BeamWidthY " << fbeamspot.BeamWidthY() << std::endl;
  }

  for (int i = 0; i<6; ++i) {
    fasciiFile << "Cov("<<i<<",j) ";
    for (int j=0; j<7; ++j) {
      fasciiFile << fbeamspot.covariance(i,j) << " ";
    }
    fasciiFile << std::endl;
  }
  // beam width error
  if (inputBeamWidth_ > 0 ) {
    fasciiFile << "Cov(6,j) 0 0 0 0 0 0 " << "1e-4" << std::endl;
  } else {
    fasciiFile << "Cov(6,j) 0 0 0 0 0 0 " << fbeamspot.covariance(6,6) << std::endl;
  }
  fasciiFile << "EmittanceX " << fbeamspot.emittanceX() << std::endl;
  fasciiFile << "EmittanceY " << fbeamspot.emittanceY() << std::endl;
  fasciiFile << "BetaStar " << fbeamspot.betaStar() << std::endl;

*/
}

void PVFitter::FitPerBunchCrossing ( )

inline

Definition at line 57 of file PVFitter.h.

References fFitPerBunchCrossing.

{ fFitPerBunchCrossing = true; }

reco::BeamSpot PVFitter::getBeamSpot ( )

inline

Definition at line 85 of file PVFitter.h.

References fbeamspot.

Referenced by AlcaBeamMonitor::endLuminosityBlock(), and BeamFitter::runPVandTrkFitter().

{ return fbeamspot; }

std::map<int, reco::BeamSpot> PVFitter::getBeamSpotMap ( )

inline

Definition at line 86 of file PVFitter.h.

References fbspotMap.

Referenced by BeamFitter::runPVandTrkFitter().

{ return fbspotMap; }

int* PVFitter::getFitLSRange ( )

inline

Definition at line 88 of file PVFitter.h.

References fbeginLumiOfFit, fendLumiOfFit, and tmp.

                       {
    int *tmp=new int[2];
    tmp[0] = fbeginLumiOfFit;
    tmp[1] = fendLumiOfFit;
    return tmp;
  }

int PVFitter::getNPVs ( )

inline

Definition at line 114 of file PVFitter.h.

References pvStore_.

Referenced by BeamFitter::getNPVs().

{ return pvStore_.size(); }

const std::map<int, int>& PVFitter::getNPVsperBX ( )

inline

Definition at line 116 of file PVFitter.h.

References bxMap_, python.multivaluedict::map(), npvsmap_, and findQualityFiles::size.

Referenced by BeamFitter::getNPVsperBX().

                                         {
    
    npvsmap_.clear();

    for ( std::map<int,std::vector<BeamSpotFitPVData> >::const_iterator pvStore = bxMap_.begin(); 
      pvStore!=bxMap_.end(); ++pvStore) {

      //std::cout << "bx " << pvStore->first << " NPVs = " << (pvStore->second).size() << std::endl;
      npvsmap_[ pvStore->first ] = (pvStore->second).size();

    }
    return npvsmap_;
  }

std::vector<BeamSpotFitPVData> PVFitter::getpvStore ( )

inline

Definition at line 55 of file PVFitter.h.

References pvStore_.

Referenced by BeamFitter::getPVvectorSize().

{ return pvStore_; } 

time_t* PVFitter::getRefTime ( )

inline

Definition at line 95 of file PVFitter.h.

References freftime.

                     {
   time_t *tmptime=new time_t[2];
   tmptime[0]=freftime[0];
   tmptime[1]=freftime[1];
   return tmptime;
  }

double PVFitter::getWidthX ( )

inline

Definition at line 48 of file PVFitter.h.

References fwidthX.

{ return fwidthX; }

double PVFitter::getWidthXerr ( )

inline

Definition at line 51 of file PVFitter.h.

References fwidthXerr.

Referenced by BeamFitter::runPVandTrkFitter().

{ return fwidthXerr; }

double PVFitter::getWidthY ( )

inline

Definition at line 49 of file PVFitter.h.

References fwidthY.

{ return fwidthY; }

double PVFitter::getWidthYerr ( )

inline

Definition at line 52 of file PVFitter.h.

References fwidthYerr.

{ return fwidthYerr; }

double PVFitter::getWidthZ ( )

inline

Definition at line 50 of file PVFitter.h.

References fwidthZ.

{ return fwidthZ; }

double PVFitter::getWidthZerr ( )

inline

Definition at line 53 of file PVFitter.h.

References fwidthZerr.

{ return fwidthZerr; }

bool PVFitter::IsFitPerBunchCrossing ( )

inline

Definition at line 87 of file PVFitter.h.

References fFitPerBunchCrossing.

Referenced by BeamFitter::dumpTxtFile(), and BeamFitter::runPVandTrkFitter().

{ return fFitPerBunchCrossing; }

double PVFitter::pvQuality

(

const reco::Vertex &

pv

)

const

vertex quality measure

Definition at line 546 of file PVFitter.cc.

References reco::Vertex::covariance().

Referenced by compressStore(), and readEvent().

{
  //
  // determinant of the transverse part of the PV covariance matrix
  //
  return
    pv.covariance(0,0)*pv.covariance(1,1)-
    pv.covariance(0,1)*pv.covariance(0,1);
}

double PVFitter::pvQuality

(

const BeamSpotFitPVData &

pv

)

const

vertex quality measure

Definition at line 557 of file PVFitter.cc.

References BeamSpotFitPVData::posCorr, and BeamSpotFitPVData::posError.

{
  //
  // determinant of the transverse part of the PV covariance matrix
  //
  double ex = pv.posError[0];
  double ey = pv.posError[1];
  return ex*ex*ey*ey*(1-pv.posCorr[0]*pv.posCorr[0]);
}

void PVFitter::readEvent

(

const edm::Event &

iEvent

)

Definition at line 83 of file PVFitter.cc.

References BeamSpotFitPVData::bunchCrossing, BeamSpotTreeData::bunchCrossing(), edm::EventBase::bunchCrossing(), bxMap_, compressStore(), dynamicQualityCut_, fFitPerBunchCrossing, ftree_, edm::Event::getByLabel(), hPVx, hPVy, edm::EventBase::id(), BeamSpotTreeData::lumi(), edm::EventBase::luminosityBlock(), maxNrVertices_, minVtxNdf_, minVtxWgt_, BeamSpotFitPVData::posCorr, BeamSpotFitPVData::posError, BeamSpotFitPVData::position, BeamSpotTreeData::pvData(), pvQuality(), pvStore_, BeamSpotTreeData::run(), edm::EventID::run(), theBeamSpotTreeData_, and vertexLabel_.

Referenced by AlcaBeamMonitor::analyze(), and BeamFitter::readEvent().

{

//   frun = iEvent.id().run();
//   const edm::TimeValue_t ftimestamp = iEvent.time().value();
//   const std::time_t ftmptime = ftimestamp >> 32;

//   if (fbeginLumiOfFit == -1) freftime[0] = freftime[1] = ftmptime;
//   if (freftime[0] == 0 || ftmptime < freftime[0]) freftime[0] = ftmptime;
//   const char* fbeginTime = formatTime(freftime[0]);
//   sprintf(fbeginTimeOfFit,"%s",fbeginTime);

//   if (freftime[1] == 0 || ftmptime > freftime[1]) freftime[1] = ftmptime;
//   const char* fendTime = formatTime(freftime[1]);
//   sprintf(fendTimeOfFit,"%s",fendTime);

//   flumi = iEvent.luminosityBlock();
//   frunFit = frun;

//   if (fbeginLumiOfFit == -1 || fbeginLumiOfFit > flumi) fbeginLumiOfFit = flumi;
//   if (fendLumiOfFit == -1 || fendLumiOfFit < flumi) fendLumiOfFit = flumi;
//   std::cout << "flumi = " <<flumi<<"; fbeginLumiOfFit = " << fbeginLumiOfFit <<"; fendLumiOfFit = "<<fendLumiOfFit<<std::endl;

  //------ Primary Vertices
  edm::Handle< reco::VertexCollection > PVCollection;
  bool hasPVs = false;
  //edm::View<reco::Vertex> vertices;
  //const reco::VertexCollection & vertices = 0;

  if ( iEvent.getByLabel(vertexLabel_, PVCollection ) ) {
      //pv = *PVCollection;
      //vertices = *PVCollection;
      hasPVs = true;
  }
  //------

  if ( hasPVs ) {

      for (reco::VertexCollection::const_iterator pv = PVCollection->begin(); pv != PVCollection->end(); ++pv ) {


           //for ( size_t ipv=0; ipv != pv.size(); ++ipv ) {

          //--- vertex selection
          if ( pv->isFake() || pv->tracksSize()==0 )  continue;
          if ( pv->ndof() < minVtxNdf_ || (pv->ndof()+3.)/pv->tracksSize()<2*minVtxWgt_ )  continue;
          //---

          hPVx->Fill( pv->x(), pv->z() );
          hPVy->Fill( pv->y(), pv->z() );

          //
          // 3D fit section
          //
          // apply additional quality cut
          if ( pvQuality(*pv)>dynamicQualityCut_ )  continue;
          // if store exceeds max. size: reduce size and apply new quality cut
          if ( pvStore_.size()>=maxNrVertices_ ) {
             compressStore();
             if ( pvQuality(*pv)>dynamicQualityCut_ )  continue;
          }
          //
          // copy PV to store
          //
      int bx = iEvent.bunchCrossing();
          BeamSpotFitPVData pvData;
      pvData.bunchCrossing = bx;
          pvData.position[0] = pv->x();
          pvData.position[1] = pv->y();
          pvData.position[2] = pv->z();
          pvData.posError[0] = pv->xError();
          pvData.posError[1] = pv->yError();
          pvData.posError[2] = pv->zError();
          pvData.posCorr[0] = pv->covariance(0,1)/pv->xError()/pv->yError();
          pvData.posCorr[1] = pv->covariance(0,2)/pv->xError()/pv->zError();
          pvData.posCorr[2] = pv->covariance(1,2)/pv->yError()/pv->zError();
          pvStore_.push_back(pvData);

      if(ftree_ != 0){
        theBeamSpotTreeData_.run(iEvent.id().run());
        theBeamSpotTreeData_.lumi(iEvent.luminosityBlock());
        theBeamSpotTreeData_.bunchCrossing(bx);
        theBeamSpotTreeData_.pvData(pvData);
        ftree_->Fill();
      }

      if (fFitPerBunchCrossing) bxMap_[bx].push_back(pvData);

      }

  }




}

void PVFitter::resetAll ( )

inline

Definition at line 74 of file PVFitter.h.

References align::BeamSpot, bxMap_, dynamicQualityCut_, fbeamspot, fbspotMap, hPVx, hPVy, pvStore_, resetLSRange(), and resetRefTime().

Referenced by BeamFitter::BeamFitter(), AlcaBeamMonitor::endLuminosityBlock(), and BeamFitter::resetPVFitter().

                  {
    resetLSRange();
    resetRefTime();
    pvStore_.clear();
    bxMap_.clear();
    dynamicQualityCut_ = 1.e30;
    hPVx->Reset();
    hPVy->Reset();
    fbeamspot = reco::BeamSpot();
    fbspotMap.clear();
  };

void PVFitter::resetLSRange ( )

inline

Definition at line 60 of file PVFitter.h.

References fbeginLumiOfFit, and fendLumiOfFit.

Referenced by resetAll().

{ fbeginLumiOfFit=fendLumiOfFit=-1; }

void PVFitter::resetRefTime ( )

inline

Definition at line 61 of file PVFitter.h.

References freftime.

Referenced by resetAll().

{ freftime[0] = freftime[1] = 0; }

void PVFitter::resizepvStore ( unsigned int  rmSize )

inline

Definition at line 103 of file PVFitter.h.

References pvStore_.

Referenced by BeamFitter::resizePVvector().

                                         {
  pvStore_.erase(pvStore_.begin(), pvStore_.begin()+rmSize);
  }

bool PVFitter::runBXFitter

(

)

Definition at line 185 of file PVFitter.cc.

References align::BeamSpot, bxMap_, errorScale_, fbeamspot, fbspotMap, fcn(), fwidthX, fwidthXerr, fwidthY, fwidthYerr, fwidthZ, fwidthZerr, python.multivaluedict::map(), makeMuonMisalignmentScenario::matrix, minNrVertices_, funct::pow(), reco::BeamSpot::setBeamWidthX(), reco::BeamSpot::setBeamWidthY(), FcnBeamSpotFitPV::setLimits(), reco::BeamSpot::setType(), sigmaCut_, findQualityFiles::size, and reco::BeamSpot::Tracker.

Referenced by BeamFitter::runPVandTrkFitter().

                           {

  edm::LogInfo("PVFitter") << " Number of bunch crossings: " << bxMap_.size() << std::endl;

  bool fit_ok = true;

  for ( std::map<int,std::vector<BeamSpotFitPVData> >::const_iterator pvStore = bxMap_.begin();
    pvStore!=bxMap_.end(); ++pvStore) {

    // first set null beam spot in case
    // fit fails
    fbspotMap[pvStore->first] = reco::BeamSpot();

    edm::LogInfo("PVFitter") << " Number of PVs collected for PVFitter: " << (pvStore->second).size() << " in bx: " << pvStore->first << std::endl;

    if ( (pvStore->second).size() <= minNrVertices_ ) {
        edm::LogWarning("PVFitter") << " not enough PVs, continue" << std::endl;
    fit_ok = false;
      continue;
    }

    //bool fit_ok = false;
    edm::LogInfo("PVFitter") << "Calculating beam spot with PVs ..." << std::endl;

    //
    // LL function and fitter
    //
    FcnBeamSpotFitPV* fcn = new FcnBeamSpotFitPV(pvStore->second);
    TFitterMinuit minuitx;
    minuitx.SetMinuitFCN(fcn);
    //
    // fit parameters: positions, widths, x-y correlations, tilts in xz and yz
    //
    minuitx.SetParameter(0,"x",0.,0.02,-10.,10.);
    minuitx.SetParameter(1,"y",0.,0.02,-10.,10.);
    minuitx.SetParameter(2,"z",0.,0.20,-30.,30.);
    minuitx.SetParameter(3,"ex",0.015,0.01,0.,10.);
    minuitx.SetParameter(4,"corrxy",0.,0.02,-1.,1.);
    minuitx.SetParameter(5,"ey",0.015,0.01,0.,10.);
    minuitx.SetParameter(6,"dxdz",0.,0.0002,-0.1,0.1);
    minuitx.SetParameter(7,"dydz",0.,0.0002,-0.1,0.1);
    minuitx.SetParameter(8,"ez",1.,0.1,0.,30.);
    minuitx.SetParameter(9,"scale",errorScale_,errorScale_/10.,errorScale_/2.,errorScale_*2.);
    //
    // first iteration without correlations
    //
    int ierr(0);
    minuitx.FixParameter(4);
    minuitx.FixParameter(6);
    minuitx.FixParameter(7);
    minuitx.FixParameter(9);
    minuitx.SetMaxIterations(100);
    //       minuitx.SetPrintLevel(3);
    minuitx.SetPrintLevel(0);
    minuitx.CreateMinimizer();
    ierr = minuitx.Minimize();
    if ( ierr ) {
        edm::LogInfo("PVFitter") << "3D beam spot fit failed in 1st iteration" << std::endl;
    fit_ok = false;
      continue;
    }
    //
    // refit with harder selection on vertices
    //
    fcn->setLimits(minuitx.GetParameter(0)-sigmaCut_*minuitx.GetParameter(3),
           minuitx.GetParameter(0)+sigmaCut_*minuitx.GetParameter(3),
           minuitx.GetParameter(1)-sigmaCut_*minuitx.GetParameter(5),
           minuitx.GetParameter(1)+sigmaCut_*minuitx.GetParameter(5),
           minuitx.GetParameter(2)-sigmaCut_*minuitx.GetParameter(8),
           minuitx.GetParameter(2)+sigmaCut_*minuitx.GetParameter(8));
    ierr = minuitx.Minimize();
    if ( ierr ) {
      edm::LogInfo("PVFitter") << "3D beam spot fit failed in 2nd iteration" << std::endl;
      fit_ok = false;
      continue;
    }
    //
    // refit with correlations
    //
    minuitx.ReleaseParameter(4);
    minuitx.ReleaseParameter(6);
    minuitx.ReleaseParameter(7);

    ierr = minuitx.Minimize();
    if ( ierr ) {
        edm::LogInfo("PVFitter") << "3D beam spot fit failed in 3rd iteration" << std::endl;
    fit_ok = false;
      continue;
    }
    // refit with floating scale factor
    //   minuitx.ReleaseParameter(9);
    //   minuitx.Minimize();

    //minuitx.PrintResults(0,0);

    fwidthX = minuitx.GetParameter(3);
    fwidthY = minuitx.GetParameter(5);
    fwidthZ = minuitx.GetParameter(8);
    fwidthXerr = minuitx.GetParError(3);
    fwidthYerr = minuitx.GetParError(5);
    fwidthZerr = minuitx.GetParError(8);

    reco::BeamSpot::CovarianceMatrix matrix;
    // need to get the full cov matrix
    matrix(0,0) = pow( minuitx.GetParError(0), 2);
    matrix(1,1) = pow( minuitx.GetParError(1), 2);
    matrix(2,2) = pow( minuitx.GetParError(2), 2);
    matrix(3,3) = fwidthZerr * fwidthZerr;
    matrix(4,4) = pow( minuitx.GetParError(6), 2);
    matrix(5,5) = pow( minuitx.GetParError(7), 2);
    matrix(6,6) = fwidthXerr * fwidthXerr;

    fbeamspot = reco::BeamSpot( reco::BeamSpot::Point(minuitx.GetParameter(0),
                              minuitx.GetParameter(1),
                              minuitx.GetParameter(2) ),
                fwidthZ,
                minuitx.GetParameter(6), minuitx.GetParameter(7),
                fwidthX,
                matrix );
    fbeamspot.setBeamWidthX( fwidthX );
    fbeamspot.setBeamWidthY( fwidthY );
    fbeamspot.setType( reco::BeamSpot::Tracker );

    fbspotMap[pvStore->first] = fbeamspot;
    edm::LogInfo("PVFitter") << "3D PV fit done for this bunch crossing."<<std::endl;
    minuitx.Clear();
    //delete fcn;
    fit_ok = fit_ok & true;
  }

  return fit_ok;
}

bool PVFitter::runFitter

(

)

Definition at line 319 of file PVFitter.cc.

References align::BeamSpot, do3DFit_, errorScale_, fbeamspot, fcn(), fwidthX, fwidthXerr, fwidthY, fwidthYerr, fwidthZ, fwidthZerr, hPVx, hPVy, edm::detail::isnan(), makeMuonMisalignmentScenario::matrix, minNrVertices_, funct::pow(), pvStore_, reco::BeamSpot::setBeamWidthX(), reco::BeamSpot::setBeamWidthY(), FcnBeamSpotFitPV::setLimits(), reco::BeamSpot::setType(), sigmaCut_, and reco::BeamSpot::Tracker.

Referenced by AlcaBeamMonitor::endLuminosityBlock(), and BeamFitter::runPVandTrkFitter().

                         {

    edm::LogInfo("PVFitter") << " Number of PVs collected for PVFitter: " << pvStore_.size() << std::endl;

    if ( pvStore_.size() <= minNrVertices_ ) return false;

    //bool fit_ok = false;

    if ( ! do3DFit_ ) {
      TH1F *h1PVx = (TH1F*) hPVx->ProjectionX("h1PVx", 0, -1, "e");
      TH1F *h1PVy = (TH1F*) hPVy->ProjectionX("h1PVy", 0, -1, "e");
      TH1F *h1PVz = (TH1F*) hPVx->ProjectionY("h1PVz", 0, -1, "e");

      h1PVx->Fit("gaus","QLM0");
      h1PVy->Fit("gaus","QLM0");
      h1PVz->Fit("gaus","QLM0");

      TF1 *gausx = h1PVx->GetFunction("gaus");
      TF1 *gausy = h1PVy->GetFunction("gaus");
      TF1 *gausz = h1PVz->GetFunction("gaus");

      fwidthX = gausx->GetParameter(2);
      fwidthY = gausy->GetParameter(2);
      fwidthZ = gausz->GetParameter(2);
      fwidthXerr = gausx->GetParError(2);
      fwidthYerr = gausy->GetParError(2);
      fwidthZerr = gausz->GetParError(2);

      reco::BeamSpot::CovarianceMatrix matrix;
      matrix(2,2) = gausz->GetParError(1) * gausz->GetParError(1);
      matrix(3,3) = fwidthZerr * fwidthZerr;
      matrix(6,6) = fwidthXerr * fwidthXerr;

      fbeamspot = reco::BeamSpot( reco::BeamSpot::Point(gausx->GetParameter(1),
                                                        gausy->GetParameter(1),
                                                        gausz->GetParameter(1) ),
                                  fwidthZ,
                                  0., 0.,
                                  fwidthX,
                                  matrix );
      fbeamspot.setBeamWidthX( fwidthX );
      fbeamspot.setBeamWidthY( fwidthY );
      fbeamspot.setType(reco::BeamSpot::Tracker);

    }
    else { // do 3D fit
      //
      // LL function and fitter
      //
      FcnBeamSpotFitPV* fcn = new FcnBeamSpotFitPV(pvStore_);
      TFitterMinuit minuitx;
      minuitx.SetMinuitFCN(fcn);
      //
      // fit parameters: positions, widths, x-y correlations, tilts in xz and yz
      //
      minuitx.SetParameter(0,"x",0.,0.02,-10.,10.);
      minuitx.SetParameter(1,"y",0.,0.02,-10.,10.);
      minuitx.SetParameter(2,"z",0.,0.20,-30.,30.);
      //minuitx.SetParameter(3,"ex",0.015,0.01,0.,10.);
      minuitx.SetParameter(3,"ex",0.015,0.01,0.0001,10.);
      minuitx.SetParameter(4,"corrxy",0.,0.02,-1.,1.);
      //minuitx.SetParameter(5,"ey",0.015,0.01,0.,10.);
      minuitx.SetParameter(5,"ey",0.015,0.01,0.0001,10.);
      minuitx.SetParameter(6,"dxdz",0.,0.0002,-0.1,0.1);
      minuitx.SetParameter(7,"dydz",0.,0.0002,-0.1,0.1);
      //minuitx.SetParameter(8,"ez",1.,0.1,0.,30.);
      minuitx.SetParameter(8,"ez",1.,0.1,1.0,30.);
      minuitx.SetParameter(9,"scale",errorScale_,errorScale_/10.,errorScale_/2.,errorScale_*2.);
      //
      // first iteration without correlations
      //
      int ierr(0);
      minuitx.FixParameter(4);
      minuitx.FixParameter(6);
      minuitx.FixParameter(7);
      minuitx.FixParameter(9);
      minuitx.SetMaxIterations(100);
//       minuitx.SetPrintLevel(3);
      minuitx.SetPrintLevel(0);
      minuitx.CreateMinimizer();
      ierr = minuitx.Minimize();
      if ( ierr ) {
          edm::LogWarning("PVFitter") << "3D beam spot fit failed in 1st iteration" << std::endl;
          return false;
      }
      //
      // refit with harder selection on vertices
      //
      fcn->setLimits(minuitx.GetParameter(0)-sigmaCut_*minuitx.GetParameter(3),
             minuitx.GetParameter(0)+sigmaCut_*minuitx.GetParameter(3),
             minuitx.GetParameter(1)-sigmaCut_*minuitx.GetParameter(5),
             minuitx.GetParameter(1)+sigmaCut_*minuitx.GetParameter(5),
             minuitx.GetParameter(2)-sigmaCut_*minuitx.GetParameter(8),
             minuitx.GetParameter(2)+sigmaCut_*minuitx.GetParameter(8));
      ierr = minuitx.Minimize();
      if ( ierr ) {
          edm::LogWarning("PVFitter") << "3D beam spot fit failed in 2nd iteration" << std::endl;
          return false;
      }
      //
      // refit with correlations
      //
      minuitx.ReleaseParameter(4);
      minuitx.ReleaseParameter(6);
      minuitx.ReleaseParameter(7);
      ierr = minuitx.Minimize();
      if ( ierr ) {
          edm::LogWarning("PVFitter") << "3D beam spot fit failed in 3rd iteration" << std::endl;
          return false;
      }
      // refit with floating scale factor
      //   minuitx.ReleaseParameter(9);
      //   minuitx.Minimize();

      //minuitx.PrintResults(0,0);

      fwidthX = minuitx.GetParameter(3);
      fwidthY = minuitx.GetParameter(5);
      fwidthZ = minuitx.GetParameter(8);
      fwidthXerr = minuitx.GetParError(3);
      fwidthYerr = minuitx.GetParError(5);
      fwidthZerr = minuitx.GetParError(8);

      // check errors on widths and sigmaZ for nan
      if ( isnan(fwidthXerr) || isnan(fwidthYerr) || isnan(fwidthZerr) ) {
          edm::LogWarning("PVFitter") << "3D beam spot fit returns nan in 3rd iteration" << std::endl;
          return false;
      }

      reco::BeamSpot::CovarianceMatrix matrix;
      // need to get the full cov matrix
      matrix(0,0) = pow( minuitx.GetParError(0), 2);
      matrix(1,1) = pow( minuitx.GetParError(1), 2);
      matrix(2,2) = pow( minuitx.GetParError(2), 2);
      matrix(3,3) = fwidthZerr * fwidthZerr;
      matrix(6,6) = fwidthXerr * fwidthXerr;

      fbeamspot = reco::BeamSpot( reco::BeamSpot::Point(minuitx.GetParameter(0),
                                                        minuitx.GetParameter(1),
                                                        minuitx.GetParameter(2) ),
                                  fwidthZ,
                                  minuitx.GetParameter(6), minuitx.GetParameter(7),
                                  fwidthX,
                                  matrix );
      fbeamspot.setBeamWidthX( fwidthX );
      fbeamspot.setBeamWidthY( fwidthY );
      fbeamspot.setType(reco::BeamSpot::Tracker);
    }

    return true; //FIXME: Need to add quality test for the fit results!
}

void PVFitter::setFitLSRange ( int  ls0, int  ls1  )

inline

Definition at line 67 of file PVFitter.h.

References fbeginLumiOfFit, and fendLumiOfFit.

                                       {
    fbeginLumiOfFit = ls0;
    fendLumiOfFit = ls1;
  }

void PVFitter::setRefTime ( std::time_t  t0, std::time_t  t1  )

inline

Definition at line 63 of file PVFitter.h.

References freftime.

                                              {
    freftime[0] = t0;
    freftime[1] = t1;
  }

void PVFitter::setTree

(

TTree *

tree

)

Definition at line 180 of file PVFitter.cc.

References BeamSpotTreeData::branch(), ftree_, theBeamSpotTreeData_, and MainPageGenerator::tree.

Referenced by BeamFitter::BeamFitter().

                                 {
  ftree_ = tree;
  theBeamSpotTreeData_.branch(ftree_);
}

Member Data Documentation

std::map< int, std::vector<BeamSpotFitPVData> > PVFitter::bxMap_

private

Definition at line 197 of file PVFitter.h.

Referenced by getNPVsperBX(), readEvent(), resetAll(), and runBXFitter().

bool PVFitter::debug_

private

Definition at line 139 of file PVFitter.h.

Referenced by PVFitter().

bool PVFitter::do3DFit_

private

Definition at line 140 of file PVFitter.h.

Referenced by PVFitter(), and runFitter().

double PVFitter::dynamicQualityCut_

private

Definition at line 198 of file PVFitter.h.

Referenced by compressStore(), PVFitter(), readEvent(), and resetAll().

double PVFitter::errorScale_

private

Definition at line 153 of file PVFitter.h.

Referenced by PVFitter(), runBXFitter(), and runFitter().

std::ofstream PVFitter::fasciiFile

private

Definition at line 137 of file PVFitter.h.

reco::BeamSpot PVFitter::fbeamspot

private

Definition at line 133 of file PVFitter.h.

Referenced by getBeamSpot(), resetAll(), runBXFitter(), and runFitter().

int PVFitter::fbeginLumiOfFit

private

Definition at line 172 of file PVFitter.h.

Referenced by getFitLSRange(), resetLSRange(), and setFitLSRange().

char PVFitter::fbeginTimeOfFit[32]

private

Definition at line 174 of file PVFitter.h.

std::map<int,reco::BeamSpot> PVFitter::fbspotMap

private

Definition at line 134 of file PVFitter.h.

Referenced by getBeamSpotMap(), resetAll(), and runBXFitter().

double PVFitter::fdxdz

private

Definition at line 187 of file PVFitter.h.

double PVFitter::fdxdzErr

private

Definition at line 193 of file PVFitter.h.

double PVFitter::fdydz

private

Definition at line 188 of file PVFitter.h.

double PVFitter::fdydzErr

private

Definition at line 194 of file PVFitter.h.

int PVFitter::fendLumiOfFit

private

Definition at line 173 of file PVFitter.h.

Referenced by getFitLSRange(), resetLSRange(), and setFitLSRange().

char PVFitter::fendTimeOfFit[32]

private

Definition at line 175 of file PVFitter.h.

bool PVFitter::fFitPerBunchCrossing

private

Definition at line 135 of file PVFitter.h.

Referenced by FitPerBunchCrossing(), IsFitPerBunchCrossing(), PVFitter(), and readEvent().

int PVFitter::flumi

private

Definition at line 157 of file PVFitter.h.

std::time_t PVFitter::freftime[2]

private

Definition at line 158 of file PVFitter.h.

Referenced by getRefTime(), resetRefTime(), and setRefTime().

int PVFitter::frun

private

Definition at line 156 of file PVFitter.h.

int PVFitter::frunFit

private

Definition at line 171 of file PVFitter.h.

double PVFitter::fsigmaZ

private

Definition at line 186 of file PVFitter.h.

double PVFitter::fsigmaZErr

private

Definition at line 192 of file PVFitter.h.

TTree* PVFitter::ftree_

private

Definition at line 162 of file PVFitter.h.

Referenced by readEvent(), and setTree().

double PVFitter::fwidthX

private

Definition at line 176 of file PVFitter.h.

Referenced by getWidthX(), runBXFitter(), and runFitter().

double PVFitter::fwidthXerr

private

Definition at line 179 of file PVFitter.h.

Referenced by getWidthXerr(), runBXFitter(), and runFitter().

double PVFitter::fwidthY

private

Definition at line 177 of file PVFitter.h.

Referenced by getWidthY(), runBXFitter(), and runFitter().

double PVFitter::fwidthYerr

private

Definition at line 180 of file PVFitter.h.

Referenced by getWidthYerr(), runBXFitter(), and runFitter().

double PVFitter::fwidthZ

private

Definition at line 178 of file PVFitter.h.

Referenced by getWidthZ(), runBXFitter(), and runFitter().

double PVFitter::fwidthZerr

private

Definition at line 181 of file PVFitter.h.

Referenced by getWidthZerr(), runBXFitter(), and runFitter().

double PVFitter::fx

private

Definition at line 183 of file PVFitter.h.

double PVFitter::fxErr

private

Definition at line 189 of file PVFitter.h.

double PVFitter::fy

private

Definition at line 184 of file PVFitter.h.

double PVFitter::fyErr

private

Definition at line 190 of file PVFitter.h.

double PVFitter::fz

private

Definition at line 185 of file PVFitter.h.

double PVFitter::fzErr

private

Definition at line 191 of file PVFitter.h.

TH2F* PVFitter::hPVx

private

Definition at line 160 of file PVFitter.h.

Referenced by PVFitter(), readEvent(), resetAll(), and runFitter().

TH2F* PVFitter::hPVy

private

Definition at line 160 of file PVFitter.h.

Referenced by PVFitter(), readEvent(), resetAll(), and runFitter().

unsigned int PVFitter::maxNrVertices_

private

Definition at line 145 of file PVFitter.h.

Referenced by PVFitter(), and readEvent().

double PVFitter::maxVtxNormChi2_

private

Definition at line 148 of file PVFitter.h.

Referenced by PVFitter().

double PVFitter::maxVtxR_

private

Definition at line 151 of file PVFitter.h.

Referenced by PVFitter().

double PVFitter::maxVtxZ_

private

Definition at line 152 of file PVFitter.h.

Referenced by PVFitter().

unsigned int PVFitter::minNrVertices_

private

Definition at line 146 of file PVFitter.h.

Referenced by PVFitter(), runBXFitter(), and runFitter().

double PVFitter::minVtxNdf_

private

Definition at line 147 of file PVFitter.h.

Referenced by PVFitter(), and readEvent().

unsigned int PVFitter::minVtxTracks_

private

Definition at line 149 of file PVFitter.h.

Referenced by PVFitter().

double PVFitter::minVtxWgt_

private

Definition at line 150 of file PVFitter.h.

Referenced by PVFitter(), and readEvent().

std::map<int, int> PVFitter::npvsmap_

private

Definition at line 132 of file PVFitter.h.

Referenced by getNPVsperBX().

std::string PVFitter::outputTxt_

private

Definition at line 143 of file PVFitter.h.

std::vector<double> PVFitter::pvQualities_

private

Definition at line 199 of file PVFitter.h.

Referenced by compressStore().

std::vector<BeamSpotFitPVData> PVFitter::pvStore_

private

Definition at line 196 of file PVFitter.h.

Referenced by compressStore(), getNPVs(), getpvStore(), readEvent(), resetAll(), resizepvStore(), and runFitter().

double PVFitter::sigmaCut_

private

Definition at line 154 of file PVFitter.h.

Referenced by PVFitter(), runBXFitter(), and runFitter().

BeamSpotTreeData PVFitter::theBeamSpotTreeData_

private

Definition at line 201 of file PVFitter.h.

Referenced by readEvent(), and setTree().

edm::InputTag PVFitter::vertexLabel_

private

Definition at line 141 of file PVFitter.h.

Referenced by PVFitter(), and readEvent().

bool PVFitter::writeTxt_

private

Definition at line 142 of file PVFitter.h.