#include <BeamFitter.h>

Public Member Functions
	BeamFitter ()
	BeamFitter (const edm::ParameterSet &iConfig)
void	dumpBWTxtFile (std::string &)
void	dumpTxtFile (std::string &, bool)
reco::BeamSpot	getBeamSpot ()
std::map< int, reco::BeamSpot >	getBeamSpotMap ()
reco::BeamSpot	getBeamWidth ()
std::vector< BSTrkParameters >	getBSvector ()
TH1F *	getCutFlow ()
int *	getFitLSRange ()
int	getNPVs ()
const std::map< int, int > &	getNPVsperBX ()
int	getNTracks ()
void	readEvent (const edm::Event &iEvent)
void	resetCutFlow ()
void	resetLSRange ()
void	resetPVFitter ()
void	resetRefTime ()
void	resetTotTrk ()
void	resetTrkVector ()
void	runAllFitter ()
bool	runBeamWidthFitter ()
bool	runFitter ()
bool	runFitterNoTxt ()
bool	runPVandTrkFitter ()
void	setFitLSRange (int ls0, int ls1)
void	setRefTime (std::time_t t0, std::time_t t1)
void	setRun (int run)
void	write2DB ()
virtual	~BeamFitter ()
Private Member Functions
const char *	formatBTime (const std::time_t &)
Private Attributes
std::vector < reco::TrackBase::TrackAlgorithm >	algorithm_
bool	appendRunTxt_
double	convergence_
int	countPass [9]
bool	debug_
bool	falgo
std::ofstream	fasciiDIP
std::ofstream	fasciiFile
reco::BeamSpot	fbeamspot
reco::BeamSpot	fbeamWidthFit
int	fbeginLumiOfFit
char	fbeginTimeOfFit [32]
std::map< int, reco::BeamSpot >	fbspotPVMap
std::vector< BSTrkParameters >	fBSvector
int	fcharge
double	fcov [7][7]
double	fd0
double	fd0bs
double	fd0phi_chi2
double	fd0phi_d0
double	fdxdz
double	fdxdzErr
double	fdydz
double	fdydzErr
int	fendLumiOfFit
char	fendTimeOfFit [32]
double	feta
bool	ffilename_changed
TFile *	file_
bool	fitted_
int	flumi
BSFitter *	fmyalgo
double	fnormchi2
int	fnPixelLayerMeas
int	fnPXBLayerMeas
int	fnPXFLayerMeas
int	fnStripLayerMeas
int	fnTECLayerMeas
int	fnTIBLayerMeas
int	fnTIDLayerMeas
int	fnTOBLayerMeas
int	fnTotLayerMeas
double	fphi0
double	fpt
TTree *	fPVTree_
bool	fpvValid
double	fpvx
double	fpvy
double	fpvz
bool	fquality
std::time_t	freftime [2]
int	frun
int	frunFit
double	fsigmad0
double	fsigmaZ
double	fsigmaz0
double	fsigmaZErr
double	ftheta
int	ftotal_tracks
TTree *	ftree_
TTree *	ftreeFit_
double	fvx
double	fvy
double	fwidthX
double	fwidthXErr
double	fwidthY
double	fwidthYErr
double	fx
double	fxErr
double	fy
double	fyErr
double	fz
double	fz0
double	fzErr
TH1F *	h1cutFlow
TH1F *	h1ntrks
TH1F *	h1vz_event
TH1F *	h1z
double	inputBeamWidth_
bool	isMuon_
int	min_Ntrks_
PVFitter *	MyPVFitter
std::string	outputDIPTxt_
std::string	outputfilename_
std::string	outputTxt_
std::vector < reco::TrackBase::TrackQuality >	quality_
bool	saveBeamFit_
bool	saveNtuple_
bool	savePVVertices_
edm::InputTag	tracksLabel_
std::vector< std::string >	trk_Algorithm_
double	trk_MaxEta_
double	trk_MaxIP_
double	trk_MaxNormChi2_
double	trk_MaxZ_
int	trk_MinNPixLayers_
int	trk_MinNTotLayers_
double	trk_MinpT_
std::vector< std::string >	trk_Quality_
edm::InputTag	vertexLabel_
bool	writeDIPTxt_
bool	writeTxt_

Detailed Description

_________________________________________________________________ class: BeamFitter.h package: RecoVertex/BeamSpotProducer

author: Francisco Yumiceva, Fermilab (yumiceva@fnal.gov) Geng-Yuan Jeng, UC Riverside (Geng-Yuan.Jeng@cern.ch)

version

Id:: BeamFitter.h,v 1.42 2010/11/08 18:53:38 friis Exp

________________________________________________________________

Definition at line 33 of file BeamFitter.h.

Constructor & Destructor Documentation

BeamFitter::BeamFitter ( ) [inline]

Definition at line 35 of file BeamFitter.h.

{}

BeamFitter::BeamFitter ( const edm::ParameterSet & iConfig )

Definition at line 42 of file BeamFitter.cc.

References reco::TrackBase::algoByName(), algorithm_, appendRunTxt_, convergence_, debug_, falgo, fasciiDIP, fbeginLumiOfFit, fbeginTimeOfFit, fBSvector, fcharge, fcov, fd0, fd0bs, fdxdz, fdxdzErr, fdydz, fdydzErr, fendLumiOfFit, fendTimeOfFit, feta, ffilename_changed, file_, fitted_, flumi, fnormchi2, fnPixelLayerMeas, fnPXBLayerMeas, fnPXFLayerMeas, fnStripLayerMeas, fnTECLayerMeas, fnTIBLayerMeas, fnTIDLayerMeas, fnTOBLayerMeas, fnTotLayerMeas, fphi0, fpt, fPVTree_, fpvValid, fpvx, fpvy, fpvz, fquality, frun, frunFit, fsigmad0, fsigmaZ, fsigmaz0, fsigmaZErr, ftheta, ftotal_tracks, ftree_, ftreeFit_, fvx, fvy, fwidthX, fwidthXErr, fwidthY, fwidthYErr, fx, fxErr, fy, fyErr, fz, fz0, fzErr, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), h1cutFlow, h1ntrks, h1vz_event, inputBeamWidth_, isMuon_, j, min_Ntrks_, MyPVFitter, outputDIPTxt_, outputfilename_, outputTxt_, ExpressReco_HICollisions_FallBack::PVFitter, quality_, reco::TrackBase::qualityByName(), PVFitter::resetAll(), resetCutFlow(), resetRefTime(), saveBeamFit_, saveNtuple_, savePVVertices_, PVFitter::setTree(), tracksLabel_, trk_Algorithm_, trk_MaxEta_, trk_MaxIP_, trk_MaxNormChi2_, trk_MaxZ_, trk_MinNPixLayers_, trk_MinNTotLayers_, trk_MinpT_, trk_Quality_, vertexLabel_, writeDIPTxt_, and writeTxt_.

                                                    : fPVTree_(0)
{

  debug_             = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("Debug");
  tracksLabel_       = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<edm::InputTag>("TrackCollection");
  vertexLabel_       = iConfig.getUntrackedParameter<edm::InputTag>("primaryVertex", edm::InputTag("offlinePrimaryVertices"));
  writeTxt_          = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("WriteAscii");
  outputTxt_         = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<std::string>("AsciiFileName");
  appendRunTxt_      = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("AppendRunToFileName");
  writeDIPTxt_       = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("WriteDIPAscii");
  outputDIPTxt_      = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<std::string>("DIPFileName");
  saveNtuple_        = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("SaveNtuple");
  saveBeamFit_       = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("SaveFitResults");
  savePVVertices_    = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("SavePVVertices");
  isMuon_            = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("IsMuonCollection");

  trk_MinpT_         = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<double>("MinimumPt");
  trk_MaxEta_        = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<double>("MaximumEta");
  trk_MaxIP_         = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<double>("MaximumImpactParameter");
  trk_MaxZ_          = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<double>("MaximumZ");
  trk_MinNTotLayers_ = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<int>("MinimumTotalLayers");
  trk_MinNPixLayers_ = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<int>("MinimumPixelLayers");
  trk_MaxNormChi2_   = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<double>("MaximumNormChi2");
  trk_Algorithm_     = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<std::vector<std::string> >("TrackAlgorithm");
  trk_Quality_       = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<std::vector<std::string> >("TrackQuality");
  min_Ntrks_         = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<int>("MinimumInputTracks");
  convergence_       = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<double>("FractionOfFittedTrks");
  inputBeamWidth_    = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<double>("InputBeamWidth",-1.);

  for (unsigned int j=0;j<trk_Algorithm_.size();j++)
    algorithm_.push_back(reco::TrackBase::algoByName(trk_Algorithm_[j]));
  for (unsigned int j=0;j<trk_Quality_.size();j++)
    quality_.push_back(reco::TrackBase::qualityByName(trk_Quality_[j]));

  if (saveNtuple_ || saveBeamFit_ || savePVVertices_) {
    outputfilename_ = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<std::string>("OutputFileName");
    file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
  }
  if (saveNtuple_) {
    ftree_ = new TTree("mytree","mytree");
    ftree_->AutoSave();

    ftree_->Branch("pt",&fpt,"fpt/D");
    ftree_->Branch("d0",&fd0,"fd0/D");
    ftree_->Branch("d0bs",&fd0bs,"fd0bs/D");
    ftree_->Branch("sigmad0",&fsigmad0,"fsigmad0/D");
    ftree_->Branch("phi0",&fphi0,"fphi0/D");
    ftree_->Branch("z0",&fz0,"fz0/D");
    ftree_->Branch("sigmaz0",&fsigmaz0,"fsigmaz0/D");
    ftree_->Branch("theta",&ftheta,"ftheta/D");
    ftree_->Branch("eta",&feta,"feta/D");
    ftree_->Branch("charge",&fcharge,"fcharge/I");
    ftree_->Branch("normchi2",&fnormchi2,"fnormchi2/D");
    ftree_->Branch("nTotLayerMeas",&fnTotLayerMeas,"fnTotLayerMeas/i");
    ftree_->Branch("nStripLayerMeas",&fnStripLayerMeas,"fnStripLayerMeas/i");
    ftree_->Branch("nPixelLayerMeas",&fnPixelLayerMeas,"fnPixelLayerMeas/i");
    ftree_->Branch("nTIBLayerMeas",&fnTIBLayerMeas,"fnTIBLayerMeas/i");
    ftree_->Branch("nTOBLayerMeas",&fnTOBLayerMeas,"fnTOBLayerMeas/i");
    ftree_->Branch("nTIDLayerMeas",&fnTIDLayerMeas,"fnTIDLayerMeas/i");
    ftree_->Branch("nTECLayerMeas",&fnTECLayerMeas,"fnTECLayerMeas/i");
    ftree_->Branch("nPXBLayerMeas",&fnPXBLayerMeas,"fnPXBLayerMeas/i");
    ftree_->Branch("nPXFLayerMeas",&fnPXFLayerMeas,"fnPXFLayerMeas/i");
    ftree_->Branch("cov",&fcov,"fcov[7][7]/D");
    ftree_->Branch("vx",&fvx,"fvx/D");
    ftree_->Branch("vy",&fvy,"fvy/D");
    ftree_->Branch("quality",&fquality,"fquality/O");
    ftree_->Branch("algo",&falgo,"falgo/O");
    ftree_->Branch("run",&frun,"frun/i");
    ftree_->Branch("lumi",&flumi,"flumi/i");
    ftree_->Branch("pvValid",&fpvValid,"fpvValid/O");
    ftree_->Branch("pvx", &fpvx, "fpvx/D");
    ftree_->Branch("pvy", &fpvy, "fpvy/D");
    ftree_->Branch("pvz", &fpvz, "fpvz/D");
  }
  if (saveBeamFit_){
    ftreeFit_ = new TTree("fitResults","fitResults");
    ftreeFit_->AutoSave();
    ftreeFit_->Branch("run",&frunFit,"frunFit/i");
    ftreeFit_->Branch("beginLumi",&fbeginLumiOfFit,"fbeginLumiOfFit/i");
    ftreeFit_->Branch("endLumi",&fendLumiOfFit,"fendLumiOfFit/i");
    ftreeFit_->Branch("beginTime",fbeginTimeOfFit,"fbeginTimeOfFit/C");
    ftreeFit_->Branch("endTime",fendTimeOfFit,"fendTimeOfFit/C");
    ftreeFit_->Branch("x",&fx,"fx/D");
    ftreeFit_->Branch("y",&fy,"fy/D");
    ftreeFit_->Branch("z",&fz,"fz/D");
    ftreeFit_->Branch("sigmaZ",&fsigmaZ,"fsigmaZ/D");
    ftreeFit_->Branch("dxdz",&fdxdz,"fdxdz/D");
    ftreeFit_->Branch("dydz",&fdydz,"fdydz/D");
    ftreeFit_->Branch("xErr",&fxErr,"fxErr/D");
    ftreeFit_->Branch("yErr",&fyErr,"fyErr/D");
    ftreeFit_->Branch("zErr",&fzErr,"fzErr/D");
    ftreeFit_->Branch("sigmaZErr",&fsigmaZErr,"fsigmaZErr/D");
    ftreeFit_->Branch("dxdzErr",&fdxdzErr,"fdxdzErr/D");
    ftreeFit_->Branch("dydzErr",&fdydzErr,"fdydzErr/D");
    ftreeFit_->Branch("widthX",&fwidthX,"fwidthX/D");
    ftreeFit_->Branch("widthY",&fwidthY,"fwidthY/D");
    ftreeFit_->Branch("widthXErr",&fwidthXErr,"fwidthXErr/D");
    ftreeFit_->Branch("widthYErr",&fwidthYErr,"fwidthYErr/D");
  }

  fBSvector.clear();
  ftotal_tracks = 0;
  fnTotLayerMeas = fnPixelLayerMeas = fnStripLayerMeas = fnTIBLayerMeas = 0;
  fnTIDLayerMeas = fnTOBLayerMeas = fnTECLayerMeas = fnPXBLayerMeas = fnPXFLayerMeas = 0;
  frun = flumi = -1;
  frunFit = fbeginLumiOfFit = fendLumiOfFit = -1;
  fquality = falgo = true;
  fpvValid = true;
  fpvx = fpvy = fpvz = 0;
  fitted_ = false;
  resetRefTime();

  //debug histograms
  h1ntrks = new TH1F("h1ntrks","number of tracks per event",50,0,50);
  h1vz_event = new TH1F("h1vz_event","track Vz", 50, -30, 30 );
  h1cutFlow = new TH1F("h1cutFlow","Cut flow table of track selection", 9, 0, 9);
  h1cutFlow->GetXaxis()->SetBinLabel(1,"No cut");
  h1cutFlow->GetXaxis()->SetBinLabel(2,"Traker hits");
  h1cutFlow->GetXaxis()->SetBinLabel(3,"Pixel hits");
  h1cutFlow->GetXaxis()->SetBinLabel(4,"norm. #chi^{2}");
  h1cutFlow->GetXaxis()->SetBinLabel(5,"algo");
  h1cutFlow->GetXaxis()->SetBinLabel(6,"quality");
  h1cutFlow->GetXaxis()->SetBinLabel(7,"d_{0}");
  h1cutFlow->GetXaxis()->SetBinLabel(8,"z_{0}");
  h1cutFlow->GetXaxis()->SetBinLabel(9,"p_{T}");
  resetCutFlow();

  // Primary vertex fitter
  MyPVFitter = new PVFitter(iConfig);
  MyPVFitter->resetAll();
  if (savePVVertices_){
    fPVTree_ = new TTree("PrimaryVertices","PrimaryVertices");
    MyPVFitter->setTree(fPVTree_);
  }

  // check filename
  ffilename_changed = false;
  if (writeDIPTxt_)
    fasciiDIP.open(outputDIPTxt_.c_str());
}

BeamFitter::~BeamFitter ( ) [virtual]

Definition at line 183 of file BeamFitter.cc.

References file_, fitted_, fPVTree_, ftree_, ftreeFit_, h1cutFlow, h1ntrks, h1vz_event, h1z, MyPVFitter, saveBeamFit_, saveNtuple_, and savePVVertices_.

                        {

  if (saveNtuple_) {
    file_->cd();
    if (fitted_ && h1z) h1z->Write();
    h1ntrks->Write();
    h1vz_event->Write();
    if (h1cutFlow) h1cutFlow->Write();
    ftree_->Write();
  }
  if (saveBeamFit_){
    file_->cd();
    ftreeFit_->Write();
  }
  if (savePVVertices_){
    file_->cd();
    fPVTree_->Write();
  }


  if (saveNtuple_ || saveBeamFit_ || savePVVertices_){
    file_->Close();
    delete file_;
  }
  delete MyPVFitter;
}

Member Function Documentation

void BeamFitter::dumpBWTxtFile ( std::string & fileName )

Definition at line 614 of file BeamFitter.cc.

References reco::BeamSpot::BeamWidthX(), reco::BeamSpot::BeamWidthXError(), fbeamWidthFit, and cmsMakeMELists::outFile.

Referenced by runBeamWidthFitter().

                                                   {
    std::ofstream outFile;
    outFile.open(fileName.c_str(),std::ios::app);
    outFile<<"-------------------------------------------------------------------------------------------------------------------------------------------------------------"<<std::endl;
    outFile<<"Beam width(in cm) from Log-likelihood fit (Here we assume a symmetric beam(SigmaX=SigmaY)!)"<<std::endl;
    outFile<<"   "<<std::endl;
    outFile << "BeamWidth =  " <<fbeamWidthFit.BeamWidthX() <<" +/- "<<fbeamWidthFit.BeamWidthXError() << std::endl;
    outFile.close();
}

void BeamFitter::dumpTxtFile	(	std::string &	fileName,
		bool	append
	)

Definition at line 624 of file BeamFitter.cc.

Referenced by runFitter(), and runPVandTrkFitter().

                                                             {
  std::ofstream outFile;

  std::string tmpname = outputTxt_;
  char index[15];
  if (appendRunTxt_ && !ffilename_changed ) {
      sprintf(index,"%s%i","_Run", frun );
      tmpname.insert(outputTxt_.length()-4,index);
      fileName = tmpname;
      ffilename_changed = true;
  }

  if (!append)
    outFile.open(fileName.c_str());
  else
    outFile.open(fileName.c_str(),std::ios::app);

  if ( MyPVFitter->IsFitPerBunchCrossing() ) {

    for (std::map<int,reco::BeamSpot>::const_iterator abspot = fbspotPVMap.begin(); abspot!= fbspotPVMap.end(); ++abspot) {
      reco::BeamSpot beamspottmp = abspot->second;
      int bx = abspot->first;

      outFile << "Runnumber " << frun << " bx " << bx << std::endl;
      outFile << "BeginTimeOfFit " << fbeginTimeOfFit << " " << freftime[0] << std::endl;
      outFile << "EndTimeOfFit " << fendTimeOfFit << " " << freftime[1] << std::endl;
      outFile << "LumiRange " << fbeginLumiOfFit << " - " << fendLumiOfFit << std::endl;
      outFile << "Type " << beamspottmp.type() << std::endl;
      outFile << "X0 " << beamspottmp.x0() << std::endl;
      outFile << "Y0 " << beamspottmp.y0() << std::endl;
      outFile << "Z0 " << beamspottmp.z0() << std::endl;
      outFile << "sigmaZ0 " << beamspottmp.sigmaZ() << std::endl;
      outFile << "dxdz " << beamspottmp.dxdz() << std::endl;
      outFile << "dydz " << beamspottmp.dydz() << std::endl;
      outFile << "BeamWidthX " << beamspottmp.BeamWidthX() << std::endl;
      outFile << "BeamWidthY " << beamspottmp.BeamWidthY() << std::endl;
      for (int i = 0; i<6; ++i) {
        outFile << "Cov("<<i<<",j) ";
        for (int j=0; j<7; ++j) {
          outFile << beamspottmp.covariance(i,j) << " ";
        }
        outFile << std::endl;
      }
      outFile << "Cov(6,j) 0 0 0 0 0 0 " << beamspottmp.covariance(6,6) << std::endl;
      //}
      outFile << "EmittanceX " << beamspottmp.emittanceX() << std::endl;
      outFile << "EmittanceY " << beamspottmp.emittanceY() << std::endl;
      outFile << "BetaStar " << beamspottmp.betaStar() << std::endl;

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

    for (int i = 0; i<6; ++i) {
      outFile << "Cov("<<i<<",j) ";
      for (int j=0; j<7; ++j) {
        outFile << fbeamspot.covariance(i,j) << " ";
      }
      outFile << std::endl;
    }

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

const char * BeamFitter::formatBTime ( const std::time_t & t ) [private]

_________________________________________________________________ class: BeamFitter.cc package: RecoVertex/BeamSpotProducer

author: Francisco Yumiceva, Fermilab (yumiceva@fnal.gov) Geng-Yuan Jeng, UC Riverside (Geng-Yuan.Jeng@cern.ch)

version

Id:: BeamFitter.cc,v 1.73 2010/11/03 15:42:09 friis Exp

________________________________________________________________

Definition at line 34 of file BeamFitter.cc.

References plotBeamSpotDB::ptm.

Referenced by runFitterNoTxt(), and runPVandTrkFitter().

                                                         {
  struct std::tm * ptm;
  ptm = gmtime(&t);
  static char ts[] = "yyyy.mn.dd hh:mm:ss zzz ";
  strftime(ts,sizeof(ts),"%Y.%m.%d %H:%M:%S %Z",ptm);
  return ts;
}

reco::BeamSpot BeamFitter::getBeamSpot ( ) [inline]

Definition at line 60 of file BeamFitter.h.

References fbeamspot.

Referenced by BeamSpotAnalyzer::endJob(), BeamSpotAnalyzer::endLuminosityBlock(), AlcaBeamSpotProducer::endLuminosityBlock(), AlcaBeamMonitor::endLuminosityBlock(), and BeamMonitor::FitAndFill().

{ return fbeamspot; }

std::map<int, reco::BeamSpot> BeamFitter::getBeamSpotMap ( ) [inline]

Definition at line 61 of file BeamFitter.h.

References fbspotPVMap.

Referenced by BeamMonitorBx::FitAndFill().

{ return fbspotPVMap; }

reco::BeamSpot BeamFitter::getBeamWidth ( ) [inline]

Definition at line 46 of file BeamFitter.h.

References fbeamWidthFit.

Referenced by BeamSpotAnalyzer::endJob().

{ return fbeamWidthFit; }

std::vector<BSTrkParameters> BeamFitter::getBSvector ( ) [inline]

Definition at line 62 of file BeamFitter.h.

References fBSvector.

Referenced by BeamSpotAnalyzer::endLuminosityBlock(), AlcaBeamSpotProducer::endLuminosityBlock(), and BeamMonitor::FitAndFill().

{ return fBSvector; }

TH1F* BeamFitter::getCutFlow ( ) [inline]

Definition at line 63 of file BeamFitter.h.

References h1cutFlow.

Referenced by BeamMonitor::analyze().

{ return h1cutFlow; }

int* BeamFitter::getFitLSRange ( ) [inline]

Definition at line 70 of file BeamFitter.h.

References fbeginLumiOfFit, fendLumiOfFit, and tmp.

Referenced by BeamSpotAnalyzer::endJob(), BeamSpotAnalyzer::endLuminosityBlock(), AlcaBeamSpotProducer::endLuminosityBlock(), BeamMonitor::FitAndFill(), and BeamMonitorBx::FitAndFill().

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

int BeamFitter::getNPVs ( ) [inline]

Definition at line 85 of file BeamFitter.h.

References PVFitter::getNPVs(), and MyPVFitter.

Referenced by BeamSpotAnalyzer::endLuminosityBlock().

                {
    return MyPVFitter->getNPVs();
  }

const std::map<int, int>& BeamFitter::getNPVsperBX ( ) [inline]

Definition at line 88 of file BeamFitter.h.

References PVFitter::getNPVsperBX(), and MyPVFitter.

Referenced by BeamMonitorBx::FitAndFill().

                                         {
    return MyPVFitter->getNPVsperBX();
  }

int BeamFitter::getNTracks ( ) [inline]

Definition at line 82 of file BeamFitter.h.

References fBSvector.

Referenced by BeamSpotAnalyzer::endLuminosityBlock().

                   {
    return fBSvector.size();
  }

void BeamFitter::readEvent ( const edm::Event & iEvent )

Definition at line 211 of file BeamFitter.cc.

Referenced by BeamMonitorBx::analyze(), BeamSpotAnalyzer::analyze(), AlcaBeamMonitor::analyze(), BeamMonitor::analyze(), and AlcaBeamSpotProducer::produce().

{


  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;
  if (freftime[1] == 0 || ftmptime > freftime[1]) freftime[1] = ftmptime;
  flumi = iEvent.luminosityBlock();
  frunFit = frun;
  if (fbeginLumiOfFit == -1 || fbeginLumiOfFit > flumi) fbeginLumiOfFit = flumi;
  if (fendLumiOfFit == -1 || fendLumiOfFit < flumi) fendLumiOfFit = flumi;

  edm::Handle<reco::TrackCollection> TrackCollection;
  iEvent.getByLabel(tracksLabel_, TrackCollection);

  //------ Primary Vertices
  edm::Handle< edm::View<reco::Vertex> > PVCollection;
  bool hasPVs = false;
  edm::View<reco::Vertex> pv;
  if ( iEvent.getByLabel(vertexLabel_, PVCollection ) ) {
      pv = *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();
  //-------

  const reco::TrackCollection *tracks = TrackCollection.product();

  double eventZ = 0;
  double averageZ = 0;

  for ( reco::TrackCollection::const_iterator track = tracks->begin();
        track != tracks->end();
        ++track ) {

    if ( ! isMuon_) {

      const reco::HitPattern& trkHP = track->hitPattern();

      fnPixelLayerMeas = trkHP.pixelLayersWithMeasurement();
      fnStripLayerMeas = trkHP.stripLayersWithMeasurement();
      fnTotLayerMeas = trkHP.trackerLayersWithMeasurement();
      fnPXBLayerMeas = trkHP.pixelBarrelLayersWithMeasurement();
      fnPXFLayerMeas = trkHP.pixelEndcapLayersWithMeasurement();
      fnTIBLayerMeas = trkHP.stripTIBLayersWithMeasurement();
      fnTIDLayerMeas = trkHP.stripTIDLayersWithMeasurement();
      fnTOBLayerMeas = trkHP.stripTOBLayersWithMeasurement();
      fnTECLayerMeas = trkHP.stripTECLayersWithMeasurement();
    } else {

      fnTotLayerMeas = track->numberOfValidHits();

    }

    fpt = track->pt();
    feta = track->eta();
    fphi0 = track->phi();
    fcharge = track->charge();
    fnormchi2 = track->normalizedChi2();
    fd0 = track->d0();
    if (refBS) fd0bs = -1*track->dxy(refBS->position());
    else fd0bs = 0.;

    fsigmad0 = track->d0Error();
    fz0 = track->dz();
    fsigmaz0 = track->dzError();
    ftheta = track->theta();
    fvx = track->vx();
    fvy = track->vy();

    for (int i=0; i<5; ++i) {
      for (int j=0; j<5; ++j) {
        fcov[i][j] = track->covariance(i,j);
      }
    }

    fquality = true;
    falgo = true;

    if (! isMuon_ ) {
      if (quality_.size()!=0) {
          fquality = false;
          for (unsigned int i = 0; i<quality_.size();++i) {
              if(debug_) edm::LogInfo("BeamFitter") << "quality_[" << i << "] = " << track->qualityName(quality_[i]) << std::endl;
              if (track->quality(quality_[i])) {
                  fquality = true;
                  break;
              }
          }
      }


      // Track algorithm

      if (algorithm_.size()!=0) {
        if (std::find(algorithm_.begin(),algorithm_.end(),track->algo())==algorithm_.end())
          falgo = false;
      }

    }

    // check if we have a valid PV
    fpvValid = false;

    if ( hasPVs ) {

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

            if (! pv[ipv].isFake() ) fpvValid = true;

            if ( ipv==0 && !pv[0].isFake() ) { fpvx = pv[0].x(); fpvy = pv[0].y(); fpvz = pv[0].z(); } // fix this later


        }

    }


    if (saveNtuple_) ftree_->Fill();
    ftotal_tracks++;

    countPass[0] = ftotal_tracks;
    // Track selection
    if (fnTotLayerMeas >= trk_MinNTotLayers_) { countPass[1] += 1;
      if (fnPixelLayerMeas >= trk_MinNPixLayers_) { countPass[2] += 1;
        if (fnormchi2 < trk_MaxNormChi2_) { countPass[3] += 1;
          if (falgo) {countPass[4] += 1;
            if (fquality) { countPass[5] += 1;
              if (std::abs( fd0 ) < trk_MaxIP_) { countPass[6] += 1;
                if (std::abs( fz0 ) < trk_MaxZ_){ countPass[7] += 1;
                  if (fpt > trk_MinpT_) {
                    countPass[8] += 1;
                    if (std::abs( feta ) < trk_MaxEta_
                        //&& fpvValid
                        ) {
                      if (debug_) {
                  edm::LogInfo("BeamFitter") << "Selected track quality = " << track->qualityMask()
                                             << "; track algorithm = " << track->algoName() << "= TrackAlgorithm: " << track->algo() << std::endl;
                      }
                      BSTrkParameters BSTrk(fz0,fsigmaz0,fd0,fsigmad0,fphi0,fpt,0.,0.);
                      BSTrk.setVx(fvx);
                      BSTrk.setVy(fvy);
                      fBSvector.push_back(BSTrk);
                      averageZ += fz0;
                    }
                  }
                }
              }
            }
          }
        }
      }
    }// track selection

  }// tracks

  averageZ = averageZ/(float)(fBSvector.size());

  for( std::vector<BSTrkParameters>::const_iterator iparam = fBSvector.begin(); iparam != fBSvector.end(); ++iparam) {

    eventZ += fabs( iparam->z0() - averageZ );

  }

  h1ntrks->Fill( fBSvector.size() );
  h1vz_event->Fill( eventZ/(float)(fBSvector.size() ) ) ;
  for (unsigned int i=0; i < sizeof(countPass)/sizeof(countPass[0]); i++)
    h1cutFlow->SetBinContent(i+1,countPass[i]);

  MyPVFitter->readEvent(iEvent);

}

void BeamFitter::resetCutFlow ( ) [inline]

Definition at line 64 of file BeamFitter.h.

References countPass, ftotal_tracks, h1cutFlow, and i.

Referenced by AlcaBeamSpotProducer::AlcaBeamSpotProducer(), BeamFitter(), BeamSpotAnalyzer::BeamSpotAnalyzer(), BeamSpotAnalyzer::endLuminosityBlock(), AlcaBeamSpotProducer::endLuminosityBlock(), BeamMonitor::FitAndFill(), and BeamMonitorBx::FitAndFill().

                      {
    h1cutFlow->Reset();
    ftotal_tracks = 0;
    for (unsigned int i=0; i<sizeof(countPass)/sizeof(countPass[0]); i++)
      countPass[i]=0;
  }

void BeamFitter::resetLSRange ( ) [inline]

Definition at line 50 of file BeamFitter.h.

References fbeginLumiOfFit, and fendLumiOfFit.

Referenced by AlcaBeamMonitor::AlcaBeamMonitor(), AlcaBeamSpotProducer::AlcaBeamSpotProducer(), BeamMonitor::BeamMonitor(), BeamMonitorBx::BeamMonitorBx(), BeamSpotAnalyzer::BeamSpotAnalyzer(), BeamSpotAnalyzer::endLuminosityBlock(), AlcaBeamSpotProducer::endLuminosityBlock(), AlcaBeamMonitor::endLuminosityBlock(), BeamMonitor::FitAndFill(), and BeamMonitorBx::FitAndFill().

{ fbeginLumiOfFit=fendLumiOfFit=-1; }

void BeamFitter::resetPVFitter ( ) [inline]

Definition at line 56 of file BeamFitter.h.

References MyPVFitter, and PVFitter::resetAll().

Referenced by AlcaBeamMonitor::AlcaBeamMonitor(), AlcaBeamSpotProducer::AlcaBeamSpotProducer(), BeamMonitor::BeamMonitor(), BeamMonitorBx::BeamMonitorBx(), BeamSpotAnalyzer::BeamSpotAnalyzer(), BeamSpotAnalyzer::endLuminosityBlock(), AlcaBeamSpotProducer::endLuminosityBlock(), AlcaBeamMonitor::endLuminosityBlock(), BeamMonitor::FitAndFill(), and BeamMonitorBx::FitAndFill().

{ MyPVFitter->resetAll(); }

void BeamFitter::resetRefTime ( ) [inline]

Definition at line 51 of file BeamFitter.h.

References freftime.

Referenced by AlcaBeamMonitor::AlcaBeamMonitor(), AlcaBeamSpotProducer::AlcaBeamSpotProducer(), BeamFitter(), BeamMonitor::BeamMonitor(), BeamMonitorBx::BeamMonitorBx(), BeamSpotAnalyzer::BeamSpotAnalyzer(), BeamSpotAnalyzer::endLuminosityBlock(), AlcaBeamSpotProducer::endLuminosityBlock(), AlcaBeamMonitor::endLuminosityBlock(), BeamMonitor::FitAndFill(), and BeamMonitorBx::FitAndFill().

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

void BeamFitter::resetTotTrk ( ) [inline]

Definition at line 49 of file BeamFitter.h.

References ftotal_tracks.

{ ftotal_tracks=0; }

void BeamFitter::resetTrkVector ( ) [inline]

Definition at line 48 of file BeamFitter.h.

References fBSvector.

Referenced by AlcaBeamMonitor::AlcaBeamMonitor(), AlcaBeamSpotProducer::AlcaBeamSpotProducer(), BeamMonitor::BeamMonitor(), BeamMonitorBx::BeamMonitorBx(), BeamSpotAnalyzer::BeamSpotAnalyzer(), BeamSpotAnalyzer::endLuminosityBlock(), AlcaBeamSpotProducer::endLuminosityBlock(), AlcaBeamMonitor::endLuminosityBlock(), BeamMonitor::FitAndFill(), and BeamMonitorBx::FitAndFill().

{ fBSvector.clear(); }

void BeamFitter::runAllFitter ( )

Definition at line 756 of file BeamFitter.cc.

References gather_cfg::cout, debug_, fbeamspot, fBSvector, BSFitter::Fit(), BSFitter::Fit_d0phi(), BSFitter::Fit_ited0phi(), BSFitter::Setd0Cut_d0phi(), BSFitter::SetFitType(), and BSFitter::SetFitVariable().

Referenced by BeamSpotAnalyzer::endJob().

                              {
  if(fBSvector.size()!=0){
    BSFitter *myalgo = new BSFitter( fBSvector );
    fbeamspot = myalgo->Fit_d0phi();

    // iterative
    if(debug_)
      std::cout << " d0-phi Iterative:" << std::endl;
    BSFitter *myitealgo = new BSFitter( fBSvector );
    myitealgo->Setd0Cut_d0phi(4.0);
    reco::BeamSpot beam_ite = myitealgo->Fit_ited0phi();
    if (debug_){
      std::cout << beam_ite << std::endl;
      std::cout << "\n Now run tests of the different fits\n";
    }
    // from here are just tests
    std::string fit_type = "chi2";
    myalgo->SetFitVariable(std::string("z"));
    myalgo->SetFitType(std::string("chi2"));
    reco::BeamSpot beam_fit_z_chi2 = myalgo->Fit();
    if (debug_){
      std::cout << " z Chi2 Fit ONLY:" << std::endl;
      std::cout << beam_fit_z_chi2 << std::endl;
    }

    fit_type = "combined";
    myalgo->SetFitVariable("z");
    myalgo->SetFitType("combined");
    reco::BeamSpot beam_fit_z_lh = myalgo->Fit();
    if (debug_){
      std::cout << " z Log-Likelihood Fit ONLY:" << std::endl;
      std::cout << beam_fit_z_lh << std::endl;
    }

    myalgo->SetFitVariable("d");
    myalgo->SetFitType("d0phi");
    reco::BeamSpot beam_fit_dphi = myalgo->Fit();
    if (debug_){
      std::cout << " d0-phi0 Fit ONLY:" << std::endl;
      std::cout << beam_fit_dphi << std::endl;
    }
    /*
    myalgo->SetFitVariable(std::string("d*z"));
    myalgo->SetFitType(std::string("likelihood"));
    reco::BeamSpot beam_fit_dz_lh = myalgo->Fit();
    if (debug_){
      std::cout << " Log-Likelihood Fit:" << std::endl;
      std::cout << beam_fit_dz_lh << std::endl;
    }

    myalgo->SetFitVariable(std::string("d*z"));
    myalgo->SetFitType(std::string("resolution"));
    reco::BeamSpot beam_fit_dresz_lh = myalgo->Fit();
    if(debug_){
      std::cout << " IP Resolution Fit" << std::endl;
      std::cout << beam_fit_dresz_lh << std::endl;

      std::cout << "c0 = " << myalgo->GetResPar0() << " +- " << myalgo->GetResPar0Err() << std::endl;
      std::cout << "c1 = " << myalgo->GetResPar1() << " +- " << myalgo->GetResPar1Err() << std::endl;
    }
    */
  }
  else
    if (debug_) std::cout << "No good track selected! No beam fit!" << std::endl;
}

bool BeamFitter::runBeamWidthFitter ( )

Definition at line 579 of file BeamFitter.cc.

References convergence_, dumpBWTxtFile(), reco::BeamSpot::Fake, fbeamspot, fbeamWidthFit, fBSvector, BSFitter::Fit(), ftotal_tracks, inputBeamWidth_, min_Ntrks_, outputTxt_, BSFitter::SetConvergence(), BSFitter::SetFitType(), BSFitter::SetFitVariable(), BSFitter::SetInputBeamWidth(), BSFitter::SetMaximumZ(), BSFitter::SetMinimumNTrks(), reco::BeamSpot::setType(), trk_MaxZ_, reco::BeamSpot::type(), and writeTxt_.

Referenced by BeamSpotAnalyzer::endJob().

                                    {
  bool widthfit_ok = false;
  // default fit to extract beam spot info
  if(fBSvector.size() > 1 ){

      edm::LogInfo("BeamFitter") << "Calculating beam spot positions('d0-phi0' method) and width using llh Fit"<< std::endl
                                 << "We will use " << fBSvector.size() << " good tracks out of " << ftotal_tracks << std::endl;

        BSFitter *myalgo = new BSFitter( fBSvector );
        myalgo->SetMaximumZ( trk_MaxZ_ );
        myalgo->SetConvergence( convergence_ );
        myalgo->SetMinimumNTrks(min_Ntrks_);
        if (inputBeamWidth_ > 0 ) myalgo->SetInputBeamWidth( inputBeamWidth_ );


   myalgo->SetFitVariable(std::string("d*z"));
   myalgo->SetFitType(std::string("likelihood"));
   fbeamWidthFit = myalgo->Fit();

   //Add to .txt file
   if(writeTxt_   ) dumpBWTxtFile(outputTxt_);

   delete myalgo;

   // not fake
   if ( fbeamspot.type() != 0 )
       widthfit_ok = true;
  }
  else{
    fbeamspot.setType(reco::BeamSpot::Fake);
    edm::LogWarning("BeamFitter") << "Not enough good tracks selected! No beam fit!" << std::endl;
  }
  return widthfit_ok;
}

bool BeamFitter::runFitter ( )

Definition at line 569 of file BeamFitter.cc.

References dumpTxtFile(), outputDIPTxt_, outputTxt_, runFitterNoTxt(), writeDIPTxt_, and writeTxt_.

                           {

    bool fit_ok = runFitterNoTxt();

    if(writeTxt_ ) dumpTxtFile(outputTxt_,true); // all reaults
    if(writeDIPTxt_) dumpTxtFile(outputDIPTxt_,false); // for DQM/DIP

    return fit_ok;
}

bool BeamFitter::runFitterNoTxt ( )

Definition at line 496 of file BeamFitter.cc.

References reco::BeamSpot::BeamWidthX(), reco::BeamSpot::BeamWidthXError(), reco::BeamSpot::BeamWidthY(), reco::BeamSpot::BeamWidthYError(), convergence_, reco::BeamSpot::dxdz(), reco::BeamSpot::dxdzError(), reco::BeamSpot::dydz(), reco::BeamSpot::dydzError(), reco::BeamSpot::Fake, fbeamspot, fbeginLumiOfFit, fbeginTimeOfFit, fBSvector, fdxdz, fdxdzErr, fdydz, fdydzErr, fendLumiOfFit, fendTimeOfFit, BSFitter::Fit(), fitted_, formatBTime(), freftime, fsigmaZ, fsigmaZErr, ftotal_tracks, ftreeFit_, fwidthX, fwidthXErr, fwidthY, fwidthYErr, fx, fxErr, fy, fyErr, fz, fzErr, BSFitter::GetVzHisto(), h1z, inputBeamWidth_, min_Ntrks_, saveBeamFit_, BSFitter::SetConvergence(), BSFitter::SetInputBeamWidth(), BSFitter::SetMaximumZ(), BSFitter::SetMinimumNTrks(), reco::BeamSpot::setType(), reco::BeamSpot::sigmaZ(), reco::BeamSpot::sigmaZ0Error(), trk_MaxZ_, reco::BeamSpot::type(), reco::BeamSpot::x0(), reco::BeamSpot::x0Error(), reco::BeamSpot::y0(), reco::BeamSpot::y0Error(), reco::BeamSpot::z0(), and reco::BeamSpot::z0Error().

Referenced by runFitter(), and runPVandTrkFitter().

                                {
  const char* fbeginTime = formatBTime(freftime[0]);
  sprintf(fbeginTimeOfFit,"%s",fbeginTime);
  const char* fendTime = formatBTime(freftime[1]);
  sprintf(fendTimeOfFit,"%s",fendTime);

  edm::LogInfo("BeamFitter") << " [BeamFitterDebugTime] freftime[0] = " << freftime[0]
                             << "; address =  " << &freftime[0]
                             << " = " << fbeginTimeOfFit << std::endl;
  edm::LogInfo("BeamFitter") << " [BeamFitterDebugTime] freftime[1] = " << freftime[1]
                             << "; address =  " << &freftime[1]
                             << " = " << fendTimeOfFit << std::endl;

  if (fbeginLumiOfFit == -1 || fendLumiOfFit == -1) {
      edm::LogWarning("BeamFitter") << "No event read! No Fitting!" << std::endl;
    return false;
  }

  bool fit_ok = false;
  // default fit to extract beam spot info
  if(fBSvector.size() > 1 ){

      edm::LogInfo("BeamFitter") << "Calculating beam spot..." << std::endl
                                 << "We will use " << fBSvector.size() << " good tracks out of " << ftotal_tracks << std::endl;

    BSFitter *myalgo = new BSFitter( fBSvector );
    myalgo->SetMaximumZ( trk_MaxZ_ );
    myalgo->SetConvergence( convergence_ );
    myalgo->SetMinimumNTrks( min_Ntrks_ );
    if (inputBeamWidth_ > 0 ) myalgo->SetInputBeamWidth( inputBeamWidth_ );

    fbeamspot = myalgo->Fit();


    // retrieve histogram for Vz
    h1z = (TH1F*) myalgo->GetVzHisto();

    delete myalgo;
    if ( fbeamspot.type() != 0 ) {// save all results except for Fake (all 0.)
      fit_ok = true;
      if (saveBeamFit_){
        fx = fbeamspot.x0();
        fy = fbeamspot.y0();
        fz = fbeamspot.z0();
        fsigmaZ = fbeamspot.sigmaZ();
        fdxdz = fbeamspot.dxdz();
        fdydz = fbeamspot.dydz();
        fwidthX = fbeamspot.BeamWidthX();
        fwidthY = fbeamspot.BeamWidthY();
        fxErr = fbeamspot.x0Error();
        fyErr = fbeamspot.y0Error();
        fzErr = fbeamspot.z0Error();
        fsigmaZErr = fbeamspot.sigmaZ0Error();
        fdxdzErr = fbeamspot.dxdzError();
        fdydzErr = fbeamspot.dydzError();
        fwidthXErr = fbeamspot.BeamWidthXError();
        fwidthYErr = fbeamspot.BeamWidthYError();
        ftreeFit_->Fill();
      }
    }
  }
  else{ // tracks <= 1
    reco::BeamSpot tmpbs;
    fbeamspot = tmpbs;
    fbeamspot.setType(reco::BeamSpot::Fake);
    edm::LogInfo("BeamFitter") << "Not enough good tracks selected! No beam fit!" << std::endl;

  }
  fitted_ = true;
  return fit_ok;

}

bool BeamFitter::runPVandTrkFitter ( )

Definition at line 394 of file BeamFitter.cc.

References reco::BeamSpot::BeamWidthX(), reco::BeamSpot::BeamWidthY(), reco::BeamSpot::covariance(), dumpTxtFile(), reco::BeamSpot::dxdz(), reco::BeamSpot::dydz(), reco::BeamSpot::Fake, fbeamspot, fbeginTimeOfFit, fbspotPVMap, fendTimeOfFit, formatBTime(), freftime, PVFitter::getBeamSpot(), PVFitter::getBeamSpotMap(), PVFitter::getWidthXerr(), inputBeamWidth_, PVFitter::IsFitPerBunchCrossing(), j, gen::k, MyPVFitter, outputDIPTxt_, outputTxt_, funct::pow(), PVFitter::runBXFitter(), PVFitter::runFitter(), runFitterNoTxt(), reco::BeamSpot::setBeamWidthY(), reco::BeamSpot::setType(), reco::BeamSpot::sigmaZ(), mathSSE::sqrt(), reco::BeamSpot::type(), reco::BeamSpot::Unknown, writeDIPTxt_, writeTxt_, reco::BeamSpot::x0(), reco::BeamSpot::y0(), and reco::BeamSpot::z0().

Referenced by BeamSpotAnalyzer::endJob(), BeamSpotAnalyzer::endLuminosityBlock(), AlcaBeamSpotProducer::endLuminosityBlock(), AlcaBeamMonitor::endLuminosityBlock(), BeamMonitor::FitAndFill(), and BeamMonitorBx::FitAndFill().

                                   {
// run both PV and track fitters
    bool fit_ok = false;
    bool pv_fit_ok = false;
    reco::BeamSpot bspotPV;
    reco::BeamSpot bspotTrk;

    // First run PV fitter
    if ( MyPVFitter->IsFitPerBunchCrossing() ){
      const char* fbeginTime = formatBTime(freftime[0]);
      sprintf(fbeginTimeOfFit,"%s",fbeginTime);
      const char* fendTime = formatBTime(freftime[1]);
      sprintf(fendTimeOfFit,"%s",fendTime);

      edm::LogInfo("BeamFitter") << " [BeamFitterBxDebugTime] freftime[0] = " << freftime[0]
                                 << "; address =  " << &freftime[0]
                                 << " = " << fbeginTimeOfFit << std::endl;
      edm::LogInfo("BeamFitter") << " [BeamFitterBxDebugTime] freftime[1] = " << freftime[1]
                                 << "; address =  " << &freftime[1]
                                 << " = " << fendTimeOfFit << std::endl;

      if ( MyPVFitter->runBXFitter() ) {
        fbspotPVMap = MyPVFitter->getBeamSpotMap();
        pv_fit_ok = true;
      }
      if(writeTxt_ ) dumpTxtFile(outputTxt_,true); // all reaults
      if(writeDIPTxt_) dumpTxtFile(outputDIPTxt_,false); // for DQM/DIP
      return pv_fit_ok;
    }

    if ( MyPVFitter->runFitter() ) {

        bspotPV = MyPVFitter->getBeamSpot();

        // take beam width from PV fit and pass it to track fitter
        // assume circular transverse beam width
        inputBeamWidth_ = sqrt( pow(bspotPV.BeamWidthX(),2) + pow(bspotPV.BeamWidthY(),2) )/sqrt(2);
        pv_fit_ok = true;

    } else {
        // problems with PV fit
        bspotPV.setType(reco::BeamSpot::Unknown);
        bspotTrk.setType(reco::BeamSpot::Unknown); //propagate error to trk beam spot
    }

    if ( runFitterNoTxt() ) {

        bspotTrk = fbeamspot;
        fit_ok = true;
    } else {
      // beamfit failed, flagged as empty beam spot
      bspotTrk.setType(reco::BeamSpot::Fake);
      fit_ok = false;
    }

    // combined results into one single beam spot

    // to pass errors I have to create another beam spot object
    reco::BeamSpot::CovarianceMatrix matrix;
    for (int j = 0 ; j < 7 ; ++j) {
        for(int k = j ; k < 7 ; ++k) {
            matrix(j,k) = bspotTrk.covariance(j,k);
        }
    }
    // change beam width error to one from PV
    if (pv_fit_ok && fit_ok ) {
      matrix(6,6) = MyPVFitter->getWidthXerr() * MyPVFitter->getWidthXerr();

      // get Z and sigmaZ from PV fit
      matrix(2,2) = bspotPV.covariance(2,2);
      matrix(3,3) = bspotPV.covariance(3,3);
      reco::BeamSpot tmpbs(reco::BeamSpot::Point(bspotTrk.x0(), bspotTrk.y0(),
                                                 bspotPV.z0() ),
                           bspotPV.sigmaZ() ,
                           bspotTrk.dxdz(),
                           bspotTrk.dydz(),
                           bspotPV.BeamWidthX(),
                           matrix,
                           bspotPV.type() );
      tmpbs.setBeamWidthY( bspotPV.BeamWidthY() );
      // overwrite beam spot result
      fbeamspot = tmpbs;
    }
    if (pv_fit_ok && fit_ok) {
      fbeamspot.setType(bspotPV.type());
    }
    else if(!pv_fit_ok && fit_ok){
      fbeamspot.setType(reco::BeamSpot::Unknown);
    }
    else if(pv_fit_ok && !fit_ok){
      fbeamspot.setType(reco::BeamSpot::Unknown);
    }
    else if(!pv_fit_ok && !fit_ok){
      fbeamspot.setType(reco::BeamSpot::Fake);
    }

    if(writeTxt_ ) dumpTxtFile(outputTxt_,true); // all reaults
    if(writeDIPTxt_) dumpTxtFile(outputDIPTxt_,false); // for DQM/DIP

    return fit_ok && pv_fit_ok;
}

void BeamFitter::setFitLSRange	(	int	ls0,
		int	ls1
	)		`[inline]`

Definition at line 76 of file BeamFitter.h.

References fbeginLumiOfFit, and fendLumiOfFit.

Referenced by BeamSpotAnalyzer::endLuminosityBlock(), AlcaBeamSpotProducer::endLuminosityBlock(), and BeamMonitor::FitAndFill().

                                      {
    fbeginLumiOfFit = ls0;
    fendLumiOfFit = ls1;
  }

void BeamFitter::setRefTime	(	std::time_t	t0,
		std::time_t	t1
	)		`[inline]`

Definition at line 52 of file BeamFitter.h.

References freftime.

Referenced by BeamSpotAnalyzer::endLuminosityBlock(), AlcaBeamSpotProducer::endLuminosityBlock(), and BeamMonitor::FitAndFill().

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

void BeamFitter::setRun ( int run ) [inline]

Definition at line 80 of file BeamFitter.h.

References frun, and CrabTask::run.

Referenced by BeamSpotAnalyzer::endLuminosityBlock(), and AlcaBeamSpotProducer::endLuminosityBlock().

{ frun = run; }

void BeamFitter::write2DB ( )

Definition at line 716 of file BeamFitter.cc.

References reco::BeamSpot::BeamWidthX(), reco::BeamSpot::BeamWidthY(), gather_cfg::cout, reco::BeamSpot::covariance(), reco::BeamSpot::dxdz(), reco::BeamSpot::dydz(), fbeamspot, i, edm::Service< T >::isAvailable(), j, reco::BeamSpot::position(), BeamSpotObjects::SetBeamWidthX(), BeamSpotObjects::SetBeamWidthY(), BeamSpotObjects::SetCovariance(), BeamSpotObjects::Setdxdz(), BeamSpotObjects::Setdydz(), BeamSpotObjects::SetPosition(), BeamSpotObjects::SetSigmaZ(), and reco::BeamSpot::sigmaZ().

Referenced by BeamSpotAnalyzer::endJob().

                         {
  BeamSpotObjects *pBSObjects = new BeamSpotObjects();

  pBSObjects->SetPosition(fbeamspot.position().X(),fbeamspot.position().Y(),fbeamspot.position().Z());
  //std::cout << " wrote: x= " << fbeamspot.position().X() << " y= "<< fbeamspot.position().Y() << " z= " << fbeamspot.position().Z() << std::endl;
  pBSObjects->SetSigmaZ(fbeamspot.sigmaZ());
  pBSObjects->Setdxdz(fbeamspot.dxdz());
  pBSObjects->Setdydz(fbeamspot.dydz());
  //if (inputBeamWidth_ > 0 ) {
  //  std::cout << " beam width value forced to be " << inputBeamWidth_ << std::endl;
  //  pBSObjects->SetBeamWidthX(inputBeamWidth_);
  //  pBSObjects->SetBeamWidthY(inputBeamWidth_);
  //} else {
    // need to fix this
    //std::cout << " using default value, 15e-4, for beam width!!!"<<std::endl;
  pBSObjects->SetBeamWidthX(fbeamspot.BeamWidthX() );
  pBSObjects->SetBeamWidthY(fbeamspot.BeamWidthY() );
  //}

  for (int i = 0; i<7; ++i) {
    for (int j=0; j<7; ++j) {
      pBSObjects->SetCovariance(i,j,fbeamspot.covariance(i,j));
    }
  }
  edm::Service<cond::service::PoolDBOutputService> poolDbService;
  if( poolDbService.isAvailable() ) {
    std::cout << "poolDBService available"<<std::endl;
    if ( poolDbService->isNewTagRequest( "BeamSpotObjectsRcd" ) ) {
      std::cout << "new tag requested" << std::endl;
      poolDbService->createNewIOV<BeamSpotObjects>( pBSObjects, poolDbService->beginOfTime(),poolDbService->endOfTime(),
                                                    "BeamSpotObjectsRcd"  );
    }
    else {
      std::cout << "no new tag requested" << std::endl;
      poolDbService->appendSinceTime<BeamSpotObjects>( pBSObjects, poolDbService->currentTime(),
                                                       "BeamSpotObjectsRcd" );
    }
  }
}