BeamFitter.cc

Go to the documentation of this file.

#include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
#include "CondFormats/BeamSpotObjects/interface/BeamSpotObjects.h"

#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Framework/interface/ConsumesCollector.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/InputTag.h"

#include "RecoVertex/BeamSpotProducer/interface/BeamFitter.h"
#include "RecoVertex/BeamSpotProducer/interface/BeamSpotWrite2Txt.h"

#include "DataFormats/Common/interface/View.h"
#include "DataFormats/TrackReco/interface/HitPattern.h"
#include "DataFormats/VertexReco/interface/Vertex.h"

// Update the string representations of the time
void BeamFitter::updateBTime() {
  char ts[] = "yyyy.mn.dd hh:mm:ss zzz ";
  char *fbeginTime = ts;
  strftime(fbeginTime, sizeof(ts), "%Y.%m.%d %H:%M:%S GMT", gmtime(&freftime[0]));
  sprintf(fbeginTimeOfFit, "%s", fbeginTime);
  char *fendTime = ts;
  strftime(fendTime, sizeof(ts), "%Y.%m.%d %H:%M:%S GMT", gmtime(&freftime[1]));
  sprintf(fendTimeOfFit, "%s", fendTime);
}

BeamFitter::BeamFitter(const edm::ParameterSet &iConfig, edm::ConsumesCollector &&iColl) : fPVTree_(nullptr) {
  debug_ = iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("Debug");
  tracksToken_ = iColl.consumes<reco::TrackCollection>(
      iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<edm::InputTag>("TrackCollection"));
  vertexToken_ = iColl.consumes<edm::View<reco::Vertex> >(
      iConfig.getUntrackedParameter<edm::InputTag>("primaryVertex", edm::InputTag("offlinePrimaryVertices")));
  beamSpotToken_ = iColl.consumes<reco::BeamSpot>(
      iConfig.getUntrackedParameter<edm::InputTag>("beamSpot", edm::InputTag("offlineBeamSpot")));
  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");
  // Specify whether we want to write the DIP file even if the fit is failed.
  writeDIPBadFit_ =
      iConfig.getParameter<edm::ParameterSet>("BeamFitter").getUntrackedParameter<bool>("WriteDIPOnBadFit", true);
  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, iColl);
  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() {
  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;
}

void BeamFitter::fillDescription(edm::ParameterSetDescription &iDesc) {
  edm::ParameterSetDescription beamFitter;

  beamFitter.addUntracked<bool>("Debug");
  beamFitter.addUntracked<edm::InputTag>("TrackCollection");
  iDesc.addUntracked<edm::InputTag>("primaryVertex", edm::InputTag("offlinePrimaryVertices"));
  iDesc.addUntracked<edm::InputTag>("beamSpot", edm::InputTag("offlineBeamSpot"));
  beamFitter.addUntracked<bool>("WriteAscii");
  beamFitter.addUntracked<std::string>("AsciiFileName");
  beamFitter.addUntracked<bool>("AppendRunToFileName");
  beamFitter.addUntracked<bool>("WriteDIPAscii");
  // Specify whether we want to write the DIP file even if the fit is failed.
  beamFitter.addUntracked<bool>("WriteDIPOnBadFit", true);
  beamFitter.addUntracked<std::string>("DIPFileName");
  beamFitter.addUntracked<bool>("SaveNtuple");
  beamFitter.addUntracked<bool>("SaveFitResults");
  beamFitter.addUntracked<bool>("SavePVVertices");
  beamFitter.addUntracked<bool>("IsMuonCollection");

  beamFitter.addUntracked<double>("MinimumPt");
  beamFitter.addUntracked<double>("MaximumEta");
  beamFitter.addUntracked<double>("MaximumImpactParameter");
  beamFitter.addUntracked<double>("MaximumZ");
  beamFitter.addUntracked<int>("MinimumTotalLayers");
  beamFitter.addUntracked<int>("MinimumPixelLayers");
  beamFitter.addUntracked<double>("MaximumNormChi2");
  beamFitter.addUntracked<std::vector<std::string> >("TrackAlgorithm");
  beamFitter.addUntracked<std::vector<std::string> >("TrackQuality");
  beamFitter.addUntracked<int>("MinimumInputTracks");
  beamFitter.addUntracked<double>("FractionOfFittedTrks");
  beamFitter.addUntracked<double>("InputBeamWidth", -1.);

  beamFitter.addUntracked<std::string>("OutputFileName", "");

  iDesc.add<edm::ParameterSetDescription>("BeamFitter", beamFitter);
}

void BeamFitter::readEvent(const edm::Event &iEvent) {
  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;
  // Update the human-readable string versions of the time
  updateBTime();

  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.getByToken(tracksToken_, TrackCollection);

  //------ Primary Vertices
  edm::Handle<edm::View<reco::Vertex> > PVCollection;
  bool hasPVs = false;
  edm::View<reco::Vertex> pv;
  if (iEvent.getByToken(vertexToken_, PVCollection)) {
    pv = *PVCollection;
    hasPVs = true;
  }
  //------

  //------ Beam Spot in current event
  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
  const reco::BeamSpot *refBS = nullptr;
  if (iEvent.getByToken(beamSpotToken_, 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_.empty()) {
        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_.empty()) {
        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);
}

bool BeamFitter::runPVandTrkFitter() {
  // 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()) {
    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_ && (pv_fit_ok || writeDIPBadFit_)) {
      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_ && ((fit_ok && pv_fit_ok) || writeDIPBadFit_)) {
    dumpTxtFile(outputDIPTxt_, false);  // for DQM/DIP
    for (size_t i = 0; i < 7; i++)
      ForDIPPV_.push_back(0.0);
  }

  return fit_ok && pv_fit_ok;
}

bool BeamFitter::runFitterNoTxt() {
  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 and Unknown (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::runFitter() {
  bool fit_ok = runFitterNoTxt();

  if (writeTxt_)
    dumpTxtFile(outputTxt_, true);  // all reaults
  if (writeDIPTxt_ && (fit_ok || writeDIPBadFit_)) {
    dumpTxtFile(outputDIPTxt_, false);  // for DQM/DIP
  }
  return fit_ok;
}

bool BeamFitter::runBeamWidthFitter() {
  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;
}

void BeamFitter::dumpBWTxtFile(std::string &fileName) {
  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) {
  std::ofstream outFile;

  std::string tmpname = outputTxt_;
  char index[15];
  if (appendRunTxt_ && writeTxt_ && !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;
    }
  }  //if bx results needed
  else {
    beamspot::BeamSpotContainer currentBS;

    currentBS.beamspot = fbeamspot;
    currentBS.run = frun;
    std::copy(fbeginTimeOfFit, fbeginTimeOfFit + 32, currentBS.beginTimeOfFit);
    std::copy(fendTimeOfFit, fendTimeOfFit + 32, currentBS.endTimeOfFit);
    currentBS.beginLumiOfFit = fbeginLumiOfFit;
    currentBS.endLumiOfFit = fendLumiOfFit;
    std::copy(freftime, freftime + 2, currentBS.reftime);

    beamspot::dumpBeamSpotTxt(outFile, currentBS);

    //write here Pv info for DIP only: This added only if append is false, which happen for DIP only :)
    if (!append) {
      outFile << "events " << (int)ForDIPPV_[0] << std::endl;
      outFile << "meanPV " << ForDIPPV_[1] << std::endl;
      outFile << "meanErrPV " << ForDIPPV_[2] << std::endl;
      outFile << "rmsPV " << ForDIPPV_[3] << std::endl;
      outFile << "rmsErrPV " << ForDIPPV_[4] << std::endl;
      outFile << "maxPV " << (int)ForDIPPV_[5] << std::endl;
      outFile << "nPV " << (int)ForDIPPV_[6] << std::endl;
    }  //writeDIPPVInfo_
  }    //else end  here

  outFile.close();
}

void BeamFitter::write2DB() {
  BeamSpotObjects pBSObjects;

  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->createOneIOV<BeamSpotObjects>(pBSObjects, poolDbService->beginOfTime(), "BeamSpotObjectsRcd");
    } else {
      std::cout << "no new tag requested" << std::endl;
      poolDbService->appendOneIOV<BeamSpotObjects>(pBSObjects, poolDbService->currentTime(), "BeamSpotObjectsRcd");
    }
  }
}

void BeamFitter::runAllFitter() {
  if (!fBSvector.empty()) {
    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;
}