CMS 3D CMS Logo

Public Member Functions | Private Attributes

PVFitter Class Reference

#include <PVFitter.h>

List of all members.

Public Member Functions

void compressStore ()
 reduce size of primary vertex cache by increasing quality limit
void dumpTxtFile ()
void FitPerBunchCrossing ()
reco::BeamSpot getBeamSpot ()
std::map< int, reco::BeamSpotgetBeamSpotMap ()
int * getFitLSRange ()
int getNPVs ()
const std::map< int, int > & getNPVsperBX ()
std::vector< BeamSpotFitPVDatagetpvStore ()
time_t * getRefTime ()
double getWidthX ()
double getWidthXerr ()
double getWidthY ()
double getWidthYerr ()
double getWidthZ ()
double getWidthZerr ()
bool IsFitPerBunchCrossing ()
 PVFitter (const edm::ParameterSet &iConfig)
 PVFitter ()
double pvQuality (const reco::Vertex &pv) const
 vertex quality measure
double pvQuality (const BeamSpotFitPVData &pv) const
 vertex quality measure
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::BeamSpotfbspotMap
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< BeamSpotFitPVDatapvStore_
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 (yumiceva@fnal.gov) Geng-Yuan Jeng, UC Riverside (Geng-Yuan.Jeng@cern.ch)

version

Id:
PVFitter.cc,v 1.18 2011/05/24 19:13:17 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 522 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 479 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.

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_, Association::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().

double PVFitter::pvQuality ( const BeamSpotFitPVData pv) const

vertex quality measure

Definition at line 565 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]);
}
double PVFitter::pvQuality ( const reco::Vertex pv) const

vertex quality measure

Definition at line 554 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);
}
void PVFitter::readEvent ( const edm::Event iEvent)

Definition at line 83 of file PVFitter.cc.

References BeamSpotTreeData::bunchCrossing(), edm::EventBase::bunchCrossing(), BeamSpotFitPVData::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, edm::Handle< T >::product(), BeamSpotTreeData::pvData(), pvQuality(), pvStore_, edm::EventID::run(), BeamSpotTreeData::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;
  }
  //------

  //------ Beam Spot in current event
  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
  const reco::BeamSpot *refBS =  0;
  if ( iEvent.getByLabel("offlineBeamSpot",recoBeamSpotHandle) )
      refBS = recoBeamSpotHandle.product();
  //-------


  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]
void PVFitter::resetLSRange ( ) [inline]

Definition at line 60 of file PVFitter.h.

References fbeginLumiOfFit, and fendLumiOfFit.

Referenced by resetAll().

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 193 of file PVFitter.cc.

References align::BeamSpot, bxMap_, errorScale_, fbeamspot, fbspotMap, fcn(), fwidthX, fwidthXerr, fwidthY, fwidthYerr, fwidthZ, fwidthZerr, Association::map, 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 327 of file PVFitter.cc.

References align::BeamSpot, do3DFit_, errorScale_, fbeamspot, fcn(), fwidthX, fwidthXerr, fwidthY, fwidthYerr, fwidthZ, fwidthZerr, hPVx, hPVy, edm::detail::isnan(), 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)

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.

Definition at line 133 of file PVFitter.h.

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

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.

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]
double PVFitter::sigmaCut_ [private]

Definition at line 154 of file PVFitter.h.

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

Definition at line 201 of file PVFitter.h.

Referenced by readEvent(), and setTree().

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.