Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-03-13 02:31:34 

0001 /*
0002  * \file BeamMonitor.cc
0003  * \author Geng-yuan Jeng/UC Riverside
0004  *         Francisco Yumiceva/FNAL
0005  */
0006 
0007 /*
0008 The code has been modified for running average
0009 mode, and it gives results for the last NLS which is
0010 configurable.
0011 Sushil S. Chauhan /UCDavis
0012 Evan Friis        /UCDavis
0013 The last tag for working versions without this change is
0014 V00-03-25
0015 */
0016 
0017 #include "DQM/BeamMonitor/plugins/BeamMonitor.h"
0018 #include "DQMServices/Core/interface/DQMStore.h"
0019 #include "FWCore/ServiceRegistry/interface/Service.h"
0020 #include "DataFormats/TrackCandidate/interface/TrackCandidate.h"
0021 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
0022 #include "DataFormats/Common/interface/View.h"
0023 #include "RecoVertex/BeamSpotProducer/interface/BSFitter.h"
0024 #include "FWCore/Framework/interface/MakerMacros.h"
0025 #include "FWCore/Framework/interface/LuminosityBlock.h"
0026 #include "FWCore/Framework/interface/Run.h"
0027 #include "FWCore/Framework/interface/ConsumesCollector.h"
0028 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0029 #include "FWCore/Common/interface/TriggerNames.h"
0030 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0031 #include "CondFormats/BeamSpotObjects/interface/BeamSpotOnlineObjects.h"
0032 #include <numeric>
0033 #include <cmath>
0034 #include <memory>
0035 #include <TMath.h>
0036 #include <iostream>
0037 #include <TStyle.h>
0038 #include <ctime>
0039 
0040 using namespace std;
0041 using namespace edm;
0042 
0043 void BeamMonitor::formatFitTime(char* ts, const time_t& t) {
0044   //constexpr int CET(+1);
0045   constexpr int CEST(+2);
0046 
0047   //tm * ptm;
0048   //ptm = gmtime ( &t );
0049   //int year = ptm->tm_year;
0050 
0051   //get correct year from ctime
0052   time_t currentTime;
0053   struct tm* localTime;
0054   time(&currentTime);                   // Get the current time
0055   localTime = localtime(&currentTime);  // Convert the current time to the local time
0056   int year = localTime->tm_year + 1900;
0057 
0058   tm* ptm;
0059   ptm = gmtime(&t);
0060 
0061   //check if year is ok
0062   if (year <= 37)
0063     year += 2000;
0064   if (year >= 70 && year <= 137)
0065     year += 1900;
0066 
0067   if (year < 1995) {
0068     edm::LogError("BadTimeStamp") << "year reported is " << year << " !!" << std::endl;
0069     //year = 2015; //overwritten later by BeamFitter.cc for fits but needed here for TH1
0070     //edm::LogError("BadTimeStamp") << "Resetting to " <<year<<std::endl;
0071   }
0072   sprintf(ts,
0073           "%4d-%02d-%02d %02d:%02d:%02d",
0074           year,
0075           ptm->tm_mon + 1,
0076           ptm->tm_mday,
0077           (ptm->tm_hour + CEST) % 24,
0078           ptm->tm_min,
0079           ptm->tm_sec);
0080 
0081 #ifdef STRIP_TRAILING_BLANKS_IN_TIMEZONE
0082   unsigned int b = strlen(ts);
0083   while (ts[--b] == ' ') {
0084     ts[b] = 0;
0085   }
0086 #endif
0087 }
0088 
0089 static constexpr int buffTime = 23;
0090 
0091 std::string BeamMonitor::getGMTstring(const time_t& timeToConvert) {
0092   char buff[32];
0093   std::strftime(buff, sizeof(buff), "%Y.%m.%d %H:%M:%S GMT", gmtime(&timeToConvert));
0094   std::string timeStr(buff);
0095   return timeStr;
0096 }
0097 
0098 //
0099 // constructors and destructor
0100 //
0101 BeamMonitor::BeamMonitor(const ParameterSet& ps)
0102     : dxBin_(ps.getParameter<int>("dxBin")),
0103       dxMin_(ps.getParameter<double>("dxMin")),
0104       dxMax_(ps.getParameter<double>("dxMax")),
0105 
0106       vxBin_(ps.getParameter<int>("vxBin")),
0107       vxMin_(ps.getParameter<double>("vxMin")),
0108       vxMax_(ps.getParameter<double>("vxMax")),
0109 
0110       phiBin_(ps.getParameter<int>("phiBin")),
0111       phiMin_(ps.getParameter<double>("phiMin")),
0112       phiMax_(ps.getParameter<double>("phiMax")),
0113 
0114       dzBin_(ps.getParameter<int>("dzBin")),
0115       dzMin_(ps.getParameter<double>("dzMin")),
0116       dzMax_(ps.getParameter<double>("dzMax")),
0117 
0118       countEvt_(0),
0119       countLumi_(0),
0120       nthBSTrk_(0),
0121       nFitElements_(3),
0122       resetHistos_(false),
0123       StartAverage_(false),
0124       firstAverageFit_(0),
0125       countGapLumi_(0),
0126       logToDb_(false),
0127       loggerActive_(false) {
0128   monitorName_ = ps.getUntrackedParameter<string>("monitorName", "YourSubsystemName");
0129   recordName_ = ps.getUntrackedParameter<string>("recordName");
0130   bsSrc_ = consumes<reco::BeamSpot>(ps.getUntrackedParameter<InputTag>("beamSpot"));
0131   tracksLabel_ = consumes<reco::TrackCollection>(
0132       ps.getParameter<ParameterSet>("BeamFitter").getUntrackedParameter<InputTag>("TrackCollection"));
0133   pvSrc_ = consumes<reco::VertexCollection>(ps.getUntrackedParameter<InputTag>("primaryVertex"));
0134   hltSrc_ = consumes<TriggerResults>(ps.getUntrackedParameter<InputTag>("hltResults"));
0135   intervalInSec_ = ps.getUntrackedParameter<int>("timeInterval", 920);  //40 LS X 23"
0136   fitNLumi_ = ps.getUntrackedParameter<int>("fitEveryNLumi", -1);
0137   resetFitNLumi_ = ps.getUntrackedParameter<int>("resetEveryNLumi", -1);
0138   fitPVNLumi_ = ps.getUntrackedParameter<int>("fitPVEveryNLumi", -1);
0139   resetPVNLumi_ = ps.getUntrackedParameter<int>("resetPVEveryNLumi", -1);
0140   deltaSigCut_ = ps.getUntrackedParameter<double>("deltaSignificanceCut", 15);
0141   debug_ = ps.getUntrackedParameter<bool>("Debug");
0142   onlineMode_ = ps.getUntrackedParameter<bool>("OnlineMode");
0143   jetTrigger_ = ps.getUntrackedParameter<std::vector<std::string> >("jetTrigger");
0144   min_Ntrks_ = ps.getParameter<ParameterSet>("BeamFitter").getUntrackedParameter<int>("MinimumInputTracks");
0145   maxZ_ = ps.getParameter<ParameterSet>("BeamFitter").getUntrackedParameter<double>("MaximumZ");
0146   minNrVertices_ = ps.getParameter<ParameterSet>("PVFitter").getUntrackedParameter<unsigned int>("minNrVerticesForFit");
0147   minVtxNdf_ = ps.getParameter<ParameterSet>("PVFitter").getUntrackedParameter<double>("minVertexNdf");
0148   minVtxWgt_ = ps.getParameter<ParameterSet>("PVFitter").getUntrackedParameter<double>("minVertexMeanWeight");
0149   useLockRecords_ = ps.getUntrackedParameter<bool>("useLockRecords");
0150   nLS_for_upload_ = ps.getUntrackedParameter<int>("nLSForUpload", 5);
0151   if (!monitorName_.empty())
0152     monitorName_ = monitorName_ + "/";
0153 
0154   theBeamFitter = std::make_unique<BeamFitter>(ps, consumesCollector());
0155   theBeamFitter->resetTrkVector();
0156   theBeamFitter->resetLSRange();
0157   theBeamFitter->resetRefTime();
0158   theBeamFitter->resetPVFitter();
0159 
0160   if (fitNLumi_ <= 0)
0161     fitNLumi_ = 1;
0162   nFits_ = beginLumiOfBSFit_ = endLumiOfBSFit_ = beginLumiOfPVFit_ = endLumiOfPVFit_ = 0;
0163   refBStime[0] = refBStime[1] = refPVtime[0] = refPVtime[1] = 0;
0164   maxZ_ = std::fabs(maxZ_);
0165   lastlumi_ = 0;
0166   nextlumi_ = 0;
0167   processed_ = false;
0168 
0169   tcdsToken_ = consumes<TCDSRecord>(ps.getUntrackedParameter<InputTag>("tcdsRecord"));
0170 }
0171 
0172 //--------------------------------------------------------
0173 namespace {
0174   /*The order of the enums is identical to the order in which
0175     MonitorElements are added to hs
0176    */
0177   enum Hists {
0178     k_x0_lumi,
0179     k_x0_lumi_all,
0180     k_y0_lumi,
0181     k_y0_lumi_all,
0182     k_z0_lumi,
0183     k_z0_lumi_all,
0184     k_sigmaX0_lumi,
0185     k_sigmaX0_lumi_all,
0186     k_sigmaY0_lumi,
0187     k_sigmaY0_lumi_all,
0188     k_sigmaZ0_lumi,
0189     k_sigmaZ0_lumi_all,
0190     k_x0_time,
0191     k_x0_time_all,
0192     k_y0_time,
0193     k_y0_time_all,
0194     k_z0_time,
0195     k_z0_time_all,
0196     k_sigmaX0_time,
0197     k_sigmaX0_time_all,
0198     k_sigmaY0_time,
0199     k_sigmaY0_time_all,
0200     k_sigmaZ0_time,
0201     k_sigmaZ0_time_all,
0202     k_PVx_lumi,
0203     k_PVx_lumi_all,
0204     k_PVy_lumi,
0205     k_PVy_lumi_all,
0206     k_PVz_lumi,
0207     k_PVz_lumi_all,
0208     k_PVx_time,
0209     k_PVx_time_all,
0210     k_PVy_time,
0211     k_PVy_time_all,
0212     k_PVz_time,
0213     k_PVz_time_all,
0214     kNumHists
0215   };
0216 }  // namespace
0217 
0218 void BeamMonitor::dqmBeginRun(edm::Run const&, edm::EventSetup const&) {
0219   if (useLockRecords_ && onlineDbService_.isAvailable()) {
0220     onlineDbService_->lockRecords();
0221   }
0222   nAnalyzedLS_ = 0;
0223 }
0224 
0225 void BeamMonitor::bookHistograms(DQMStore::IBooker& iBooker, edm::Run const& iRun, edm::EventSetup const& iSetup) {
0226   frun = iRun.run();
0227   ftimestamp = iRun.beginTime().value();
0228   tmpTime = ftimestamp >> 32;
0229   startTime = refTime = tmpTime;
0230   char eventTime[64];
0231   formatFitTime(eventTime, tmpTime);
0232   edm::LogInfo("BeamMonitor") << "TimeOffset = " << eventTime << std::endl;
0233   TDatime da(eventTime);
0234   if (debug_) {
0235     edm::LogInfo("BeamMonitor") << "TimeOffset = ";
0236     da.Print();
0237   }
0238   auto daTime = da.Convert(kTRUE);
0239 
0240   // book some histograms here
0241 
0242   // create and cd into new folder
0243   iBooker.setCurrentFolder(monitorName_ + "Fit");
0244 
0245   h_nTrk_lumi = iBooker.book1D("nTrk_lumi", "Num. of selected tracks vs lumi (Fit)", 20, 0.5, 20.5);
0246   h_nTrk_lumi->setAxisTitle("Lumisection", 1);
0247   h_nTrk_lumi->setAxisTitle("Num of Tracks for Fit", 2);
0248 
0249   //store vtx vs lumi for monitoring why fits fail
0250   h_nVtx_lumi = iBooker.book1D("nVtx_lumi", "Num. of selected Vtx vs lumi (Fit)", 20, 0.5, 20.5);
0251   h_nVtx_lumi->setAxisTitle("Lumisection", 1);
0252   h_nVtx_lumi->setAxisTitle("Num of Vtx for Fit", 2);
0253 
0254   h_nVtx_lumi_all = iBooker.book1D("nVtx_lumi_all", "Num. of selected Vtx vs lumi (Fit) all", 20, 0.5, 20.5);
0255   h_nVtx_lumi_all->getTH1()->SetCanExtend(TH1::kAllAxes);
0256   h_nVtx_lumi_all->setAxisTitle("Lumisection", 1);
0257   h_nVtx_lumi_all->setAxisTitle("Num of Vtx for Fit", 2);
0258 
0259   h_d0_phi0 = iBooker.bookProfile(
0260       "d0_phi0", "d_{0} vs. #phi_{0} (Selected Tracks)", phiBin_, phiMin_, phiMax_, dxBin_, dxMin_, dxMax_, "");
0261   h_d0_phi0->setAxisTitle("#phi_{0} (rad)", 1);
0262   h_d0_phi0->setAxisTitle("d_{0} (cm)", 2);
0263 
0264   h_vx_vy = iBooker.book2D(
0265       "trk_vx_vy", "Vertex (PCA) position of selected tracks", vxBin_, vxMin_, vxMax_, vxBin_, vxMin_, vxMax_);
0266   h_vx_vy->setOption("COLZ");
0267   //   h_vx_vy->getTH1()->SetCanExtend(TH1::kAllAxes);
0268   h_vx_vy->setAxisTitle("x coordinate of input track at PCA (cm)", 1);
0269   h_vx_vy->setAxisTitle("y coordinate of input track at PCA (cm)", 2);
0270 
0271   {
0272     TDatime* da = new TDatime();
0273     gStyle->SetTimeOffset(da->Convert(kTRUE));
0274   }
0275 
0276   const int nvar_ = 6;
0277   string coord[nvar_] = {"x", "y", "z", "sigmaX", "sigmaY", "sigmaZ"};
0278   string label[nvar_] = {
0279       "x_{0} (cm)", "y_{0} (cm)", "z_{0} (cm)", "#sigma_{X_{0}} (cm)", "#sigma_{Y_{0}} (cm)", "#sigma_{Z_{0}} (cm)"};
0280 
0281   hs.reserve(kNumHists);
0282   for (int i = 0; i < 4; i++) {
0283     iBooker.setCurrentFolder(monitorName_ + "Fit");
0284     for (int ic = 0; ic < nvar_; ++ic) {
0285       TString histName(coord[ic]);
0286       TString histTitle(coord[ic]);
0287       string ytitle(label[ic]);
0288       string xtitle("");
0289       string options("E1");
0290       bool createHisto = true;
0291       switch (i) {
0292         case 1:  // BS vs time
0293           histName += "0_time";
0294           xtitle = "Time [UTC]";
0295           if (ic < 3)
0296             histTitle += " coordinate of beam spot vs time (Fit)";
0297           else
0298             histTitle = histTitle.Insert(5, " ") + " of beam spot vs time (Fit)";
0299           break;
0300         case 2:  // PV vs lumi
0301           if (ic < 3) {
0302             iBooker.setCurrentFolder(monitorName_ + "PrimaryVertex");
0303             histName.Insert(0, "PV");
0304             histName += "_lumi";
0305             histTitle.Insert(0, "Avg. ");
0306             histTitle += " position of primary vtx vs lumi";
0307             xtitle = "Lumisection";
0308             ytitle.insert(0, "PV");
0309             ytitle += " #pm #sigma_{PV";
0310             ytitle += coord[ic];
0311             ytitle += "} (cm)";
0312           } else
0313             createHisto = false;
0314           break;
0315         case 3:  // PV vs time
0316           if (ic < 3) {
0317             iBooker.setCurrentFolder(monitorName_ + "PrimaryVertex");
0318             histName.Insert(0, "PV");
0319             histName += "_time";
0320             histTitle.Insert(0, "Avg. ");
0321             histTitle += " position of primary vtx vs time";
0322             xtitle = "Time [UTC]";
0323             ytitle.insert(0, "PV");
0324             ytitle += " #pm #sigma_{PV";
0325             ytitle += coord[ic];
0326             ytitle += "} (cm)";
0327           } else
0328             createHisto = false;
0329           break;
0330         default:  // BS vs lumi
0331           histName += "0_lumi";
0332           xtitle = "Lumisection";
0333           if (ic < 3)
0334             histTitle += " coordinate of beam spot vs lumi (Fit)";
0335           else
0336             histTitle = histTitle.Insert(5, " ") + " of beam spot vs lumi (Fit)";
0337           break;
0338       }
0339       if (createHisto) {
0340         edm::LogInfo("BeamMonitor") << "hitsName = " << histName << "; histTitle = " << histTitle << std::endl;
0341         auto tmpHs = iBooker.book1D(histName, histTitle, 40, 0.5, 40.5);
0342         hs.push_back(tmpHs);
0343         tmpHs->setAxisTitle(xtitle, 1);
0344         tmpHs->setAxisTitle(ytitle, 2);
0345         tmpHs->getTH1()->SetOption("E1");
0346         if (histName.Contains("time")) {
0347           //int nbins = (intervalInSec_/23 > 0 ? intervalInSec_/23 : 40);
0348           tmpHs->getTH1()->SetBins(intervalInSec_, 0.5, intervalInSec_ + 0.5);
0349           tmpHs->setAxisTimeDisplay(1);
0350           tmpHs->setAxisTimeFormat("%H:%M:%S", 1);
0351           tmpHs->getTH1()->GetXaxis()->SetTimeOffset(daTime);
0352         }
0353         histName += "_all";
0354         histTitle += " all";
0355         tmpHs = iBooker.book1D(histName, histTitle, 40, 0.5, 40.5);
0356         hs.push_back(tmpHs);
0357         tmpHs->getTH1()->SetCanExtend(TH1::kAllAxes);
0358         tmpHs->setAxisTitle(xtitle, 1);
0359         tmpHs->setAxisTitle(ytitle, 2);
0360         tmpHs->getTH1()->SetOption("E1");
0361         if (histName.Contains("time")) {
0362           //int nbins = (intervalInSec_/23 > 0 ? intervalInSec_/23 : 40);
0363           tmpHs->getTH1()->SetBins(intervalInSec_, 0.5, intervalInSec_ + 0.5);
0364           tmpHs->setAxisTimeDisplay(1);
0365           tmpHs->setAxisTimeFormat("%H:%M:%S", 1);
0366           tmpHs->getTH1()->GetXaxis()->SetTimeOffset(daTime);
0367         }
0368       }
0369     }
0370   }
0371   assert(hs.size() == kNumHists);
0372   assert(0 == strcmp(hs[k_sigmaY0_time]->getTH1()->GetName(), "sigmaY0_time"));
0373   assert(0 == strcmp(hs[k_PVz_lumi_all]->getTH1()->GetName(), "PVz_lumi_all"));
0374 
0375   iBooker.setCurrentFolder(monitorName_ + "Fit");
0376 
0377   h_trk_z0 = iBooker.book1D("trk_z0", "z_{0} of selected tracks", dzBin_, dzMin_, dzMax_);
0378   h_trk_z0->setAxisTitle("z_{0} of selected tracks (cm)", 1);
0379 
0380   h_vx_dz = iBooker.bookProfile(
0381       "vx_dz", "v_{x} vs. dz of selected tracks", dzBin_, dzMin_, dzMax_, dxBin_, dxMin_, dxMax_, "");
0382   h_vx_dz->setAxisTitle("dz (cm)", 1);
0383   h_vx_dz->setAxisTitle("x coordinate of input track at PCA (cm)", 2);
0384 
0385   h_vy_dz = iBooker.bookProfile(
0386       "vy_dz", "v_{y} vs. dz of selected tracks", dzBin_, dzMin_, dzMax_, dxBin_, dxMin_, dxMax_, "");
0387   h_vy_dz->setAxisTitle("dz (cm)", 1);
0388   h_vy_dz->setAxisTitle("y coordinate of input track at PCA (cm)", 2);
0389 
0390   h_x0 = iBooker.book1D("BeamMonitorFeedBack_x0", "x coordinate of beam spot (Fit)", 100, -0.01, 0.01);
0391   h_x0->setAxisTitle("x_{0} (cm)", 1);
0392   h_x0->getTH1()->SetCanExtend(TH1::kAllAxes);
0393 
0394   h_y0 = iBooker.book1D("BeamMonitorFeedBack_y0", "y coordinate of beam spot (Fit)", 100, -0.01, 0.01);
0395   h_y0->setAxisTitle("y_{0} (cm)", 1);
0396   h_y0->getTH1()->SetCanExtend(TH1::kAllAxes);
0397 
0398   h_z0 = iBooker.book1D("BeamMonitorFeedBack_z0", "z coordinate of beam spot (Fit)", dzBin_, dzMin_, dzMax_);
0399   h_z0->setAxisTitle("z_{0} (cm)", 1);
0400   h_z0->getTH1()->SetCanExtend(TH1::kAllAxes);
0401 
0402   h_sigmaX0 = iBooker.book1D("BeamMonitorFeedBack_sigmaX0", "sigma x0 of beam spot (Fit)", 100, 0, 0.05);
0403   h_sigmaX0->setAxisTitle("#sigma_{X_{0}} (cm)", 1);
0404   h_sigmaX0->getTH1()->SetCanExtend(TH1::kAllAxes);
0405 
0406   h_sigmaY0 = iBooker.book1D("BeamMonitorFeedBack_sigmaY0", "sigma y0 of beam spot (Fit)", 100, 0, 0.05);
0407   h_sigmaY0->setAxisTitle("#sigma_{Y_{0}} (cm)", 1);
0408   h_sigmaY0->getTH1()->SetCanExtend(TH1::kAllAxes);
0409 
0410   h_sigmaZ0 = iBooker.book1D("BeamMonitorFeedBack_sigmaZ0", "sigma z0 of beam spot (Fit)", 100, 0, 10);
0411   h_sigmaZ0->setAxisTitle("#sigma_{Z_{0}} (cm)", 1);
0412   h_sigmaZ0->getTH1()->SetCanExtend(TH1::kAllAxes);
0413 
0414   // Histograms of all reco tracks (without cuts):
0415   h_trkPt = iBooker.book1D("trkPt", "p_{T} of all reco'd tracks (no selection)", 200, 0., 50.);
0416   h_trkPt->setAxisTitle("p_{T} (GeV/c)", 1);
0417 
0418   h_trkVz = iBooker.book1D("trkVz", "Z coordinate of PCA of all reco'd tracks (no selection)", dzBin_, dzMin_, dzMax_);
0419   h_trkVz->setAxisTitle("V_{Z} (cm)", 1);
0420 
0421   cutFlowTable = iBooker.book1D("cutFlowTable", "Cut flow table of track selection", 9, 0, 9);
0422 
0423   // Results of previous good fit:
0424   fitResults = iBooker.book2D("fitResults", "Results of previous good beam fit", 2, 0, 2, 8, 0, 8);
0425   fitResults->setAxisTitle("Fitted Beam Spot (cm)", 1);
0426   fitResults->setBinLabel(8, "x_{0}", 2);
0427   fitResults->setBinLabel(7, "y_{0}", 2);
0428   fitResults->setBinLabel(6, "z_{0}", 2);
0429   fitResults->setBinLabel(5, "#sigma_{Z}", 2);
0430   fitResults->setBinLabel(4, "#frac{dx}{dz} (rad)", 2);
0431   fitResults->setBinLabel(3, "#frac{dy}{dz} (rad)", 2);
0432   fitResults->setBinLabel(2, "#sigma_{X}", 2);
0433   fitResults->setBinLabel(1, "#sigma_{Y}", 2);
0434   fitResults->setBinLabel(1, "Mean", 1);
0435   fitResults->setBinLabel(2, "Stat. Error", 1);
0436   fitResults->getTH1()->SetOption("text");
0437 
0438   // Histos of PrimaryVertices:
0439   iBooker.setCurrentFolder(monitorName_ + "PrimaryVertex");
0440 
0441   h_nVtx = iBooker.book1D("vtxNbr", "Reconstructed Vertices(non-fake) in all Event", 60, -0.5, 59.5);
0442   h_nVtx->setAxisTitle("Num. of reco. vertices", 1);
0443 
0444   //For one Trigger only
0445   h_nVtx_st = iBooker.book1D("vtxNbr_SelectedTriggers", "Reconstructed Vertices(non-fake) in Events", 60, -0.5, 59.5);
0446   //h_nVtx_st->setAxisTitle("Num. of reco. vertices for Un-Prescaled Jet Trigger",1);
0447 
0448   // Monitor only the PV with highest sum pt of assoc. trks:
0449   h_PVx[0] = iBooker.book1D("PVX", "x coordinate of Primary Vtx", 50, -0.01, 0.01);
0450   h_PVx[0]->setAxisTitle("PVx (cm)", 1);
0451   h_PVx[0]->getTH1()->SetCanExtend(TH1::kAllAxes);
0452 
0453   h_PVy[0] = iBooker.book1D("PVY", "y coordinate of Primary Vtx", 50, -0.01, 0.01);
0454   h_PVy[0]->setAxisTitle("PVy (cm)", 1);
0455   h_PVy[0]->getTH1()->SetCanExtend(TH1::kAllAxes);
0456 
0457   h_PVz[0] = iBooker.book1D("PVZ", "z coordinate of Primary Vtx", dzBin_, dzMin_, dzMax_);
0458   h_PVz[0]->setAxisTitle("PVz (cm)", 1);
0459 
0460   h_PVx[1] = iBooker.book1D("PVXFit", "x coordinate of Primary Vtx (Last Fit)", 50, -0.01, 0.01);
0461   h_PVx[1]->setAxisTitle("PVx (cm)", 1);
0462   h_PVx[1]->getTH1()->SetCanExtend(TH1::kAllAxes);
0463 
0464   h_PVy[1] = iBooker.book1D("PVYFit", "y coordinate of Primary Vtx (Last Fit)", 50, -0.01, 0.01);
0465   h_PVy[1]->setAxisTitle("PVy (cm)", 1);
0466   h_PVy[1]->getTH1()->SetCanExtend(TH1::kAllAxes);
0467 
0468   h_PVz[1] = iBooker.book1D("PVZFit", "z coordinate of Primary Vtx (Last Fit)", dzBin_, dzMin_, dzMax_);
0469   h_PVz[1]->setAxisTitle("PVz (cm)", 1);
0470 
0471   h_PVxz = iBooker.bookProfile("PVxz", "PVx vs. PVz", dzBin_ / 2, dzMin_, dzMax_, dxBin_ / 2, dxMin_, dxMax_, "");
0472   h_PVxz->setAxisTitle("PVz (cm)", 1);
0473   h_PVxz->setAxisTitle("PVx (cm)", 2);
0474 
0475   h_PVyz = iBooker.bookProfile("PVyz", "PVy vs. PVz", dzBin_ / 2, dzMin_, dzMax_, dxBin_ / 2, dxMin_, dxMax_, "");
0476   h_PVyz->setAxisTitle("PVz (cm)", 1);
0477   h_PVyz->setAxisTitle("PVy (cm)", 2);
0478 
0479   // Results of previous good fit:
0480   pvResults = iBooker.book2D("pvResults", "Results of fitting Primary Vertices", 2, 0, 2, 6, 0, 6);
0481   pvResults->setAxisTitle("Fitted Primary Vertex (cm)", 1);
0482   pvResults->setBinLabel(6, "PVx", 2);
0483   pvResults->setBinLabel(5, "PVy", 2);
0484   pvResults->setBinLabel(4, "PVz", 2);
0485   pvResults->setBinLabel(3, "#sigma_{X}", 2);
0486   pvResults->setBinLabel(2, "#sigma_{Y}", 2);
0487   pvResults->setBinLabel(1, "#sigma_{Z}", 2);
0488   pvResults->setBinLabel(1, "Mean", 1);
0489   pvResults->setBinLabel(2, "Stat. Error", 1);
0490   pvResults->getTH1()->SetOption("text");
0491 
0492   // Summary plots:
0493   iBooker.setCurrentFolder(monitorName_ + "EventInfo");
0494 
0495   reportSummary = iBooker.bookFloat("reportSummary");
0496   if (reportSummary)
0497     reportSummary->Fill(std::numeric_limits<double>::quiet_NaN());
0498 
0499   char histo[20];
0500   iBooker.setCurrentFolder(monitorName_ + "EventInfo/reportSummaryContents");
0501   for (int n = 0; n < nFitElements_; n++) {
0502     switch (n) {
0503       case 0:
0504         sprintf(histo, "x0_status");
0505         break;
0506       case 1:
0507         sprintf(histo, "y0_status");
0508         break;
0509       case 2:
0510         sprintf(histo, "z0_status");
0511         break;
0512     }
0513     reportSummaryContents[n] = iBooker.bookFloat(histo);
0514   }
0515 
0516   for (int i = 0; i < nFitElements_; i++) {
0517     summaryContent_[i] = 0.;
0518     reportSummaryContents[i]->Fill(std::numeric_limits<double>::quiet_NaN());
0519   }
0520 
0521   iBooker.setCurrentFolder(monitorName_ + "EventInfo");
0522 
0523   reportSummaryMap = iBooker.book2D("reportSummaryMap", "Beam Spot Summary Map", 1, 0, 1, 3, 0, 3);
0524   reportSummaryMap->setAxisTitle("", 1);
0525   reportSummaryMap->setAxisTitle("Fitted Beam Spot", 2);
0526   reportSummaryMap->setBinLabel(1, " ", 1);
0527   reportSummaryMap->setBinLabel(1, "x_{0}", 2);
0528   reportSummaryMap->setBinLabel(2, "y_{0}", 2);
0529   reportSummaryMap->setBinLabel(3, "z_{0}", 2);
0530   for (int i = 0; i < nFitElements_; i++) {
0531     reportSummaryMap->setBinContent(1, i + 1, -1.);
0532   }
0533 }
0534 
0535 //--------------------------------------------------------
0536 void BeamMonitor::beginLuminosityBlock(const LuminosityBlock& lumiSeg, const EventSetup& context) {
0537   // start DB logger
0538   DBloggerReturn_ = 0;
0539   nAnalyzedLS_++;
0540   if (onlineDbService_.isAvailable() && logToDb_) {
0541     onlineDbService_->logger().start();
0542     loggerActive_ = true;
0543     onlineDbService_->logger().logInfo() << "BeamMonitor::beginLuminosityBlock - LS: " << lumiSeg.luminosityBlock()
0544                                          << " - Run: " << lumiSeg.getRun().run();
0545   }
0546 
0547   int nthlumi = lumiSeg.luminosityBlock();
0548   const edm::TimeValue_t fbegintimestamp = lumiSeg.beginTime().value();
0549   const std::time_t ftmptime = fbegintimestamp >> 32;
0550 
0551   if (countLumi_ == 0 && (!processed_)) {
0552     beginLumiOfBSFit_ = beginLumiOfPVFit_ = nthlumi;
0553     refBStime[0] = refPVtime[0] = ftmptime;
0554     mapBeginBSLS[countLumi_] = nthlumi;
0555     mapBeginPVLS[countLumi_] = nthlumi;
0556     mapBeginBSTime[countLumi_] = ftmptime;
0557     mapBeginPVTime[countLumi_] = ftmptime;
0558   }  //for the first record
0559 
0560   if (nthlumi > nextlumi_) {
0561     if (processed_) {
0562       countLumi_++;
0563       //store here them will need when we remove the first one of Last N LS
0564       mapBeginBSLS[countLumi_] = nthlumi;
0565       mapBeginPVLS[countLumi_] = nthlumi;
0566       mapBeginBSTime[countLumi_] = ftmptime;
0567       mapBeginPVTime[countLumi_] = ftmptime;
0568     }  //processed passed but not the first lumi
0569     if ((!processed_) && countLumi_ != 0) {
0570       mapBeginBSLS[countLumi_] = nthlumi;
0571       mapBeginPVLS[countLumi_] = nthlumi;
0572       mapBeginBSTime[countLumi_] = ftmptime;
0573       mapBeginPVTime[countLumi_] = ftmptime;
0574     }  //processed fails for last lumi
0575   }  //nthLumi > nextlumi
0576 
0577   if (StartAverage_) {
0578     //Just Make sure it get rest
0579     refBStime[0] = 0;
0580     refPVtime[0] = 0;
0581     beginLumiOfPVFit_ = 0;
0582     beginLumiOfBSFit_ = 0;
0583 
0584     if (debug_)
0585       edm::LogInfo("BeamMonitor") << " beginLuminosityBlock:  Size of mapBeginBSLS before =  " << mapBeginBSLS.size()
0586                                   << endl;
0587     if (nthlumi >
0588         nextlumi_) {  //this make sure that it does not take into account this lumi for fitting and only look forward for new lumi
0589       //as countLumi also remains the same so map value  get overwritten once return to normal running.
0590       //even if few LS are misssing and DQM module do not sees them then it catchs up again
0591       map<int, int>::iterator itbs = mapBeginBSLS.begin();
0592       map<int, int>::iterator itpv = mapBeginPVLS.begin();
0593       map<int, std::time_t>::iterator itbstime = mapBeginBSTime.begin();
0594       map<int, std::time_t>::iterator itpvtime = mapBeginPVTime.begin();
0595 
0596       if (processed_) {  // otherwise if false then LS range of fit get messed up because we don't remove trk/pvs but we remove LS begin value . This prevent it as it happened if LS is there but no event are processed for some reason
0597         mapBeginBSLS.erase(itbs);
0598         mapBeginPVLS.erase(itpv);
0599         mapBeginBSTime.erase(itbstime);
0600         mapBeginPVTime.erase(itpvtime);
0601       }
0602       /*//not sure if want this or not ??
0603             map<int, int>::iterator itgapb=mapBeginBSLS.begin();
0604             map<int, int>::iterator itgape=mapBeginBSLS.end(); itgape--;
0605             countGapLumi_ = ( (itgape->second) - (itgapb->second) );
0606             //if we see Gap more than then 2*resetNFitLumi !!!!!!!
0607             //for example if 10-15 is fitted and if 16-25 are missing then we next fit will be for range 11-26 but BS can change in between
0608             // so better start  as fresh  and reset everything like starting in the begining!
0609             if(countGapLumi_ >= 2*resetFitNLumi_){RestartFitting(); mapBeginBSLS[countLumi_]   = nthlumi;}
0610             */
0611     }
0612 
0613     if (debug_)
0614       edm::LogInfo("BeamMonitor") << " beginLuminosityBlock::  Size of mapBeginBSLS After = " << mapBeginBSLS.size()
0615                                   << endl;
0616 
0617     map<int, int>::iterator bbs = mapBeginBSLS.begin();
0618     map<int, int>::iterator bpv = mapBeginPVLS.begin();
0619     map<int, std::time_t>::iterator bbst = mapBeginBSTime.begin();
0620     map<int, std::time_t>::iterator bpvt = mapBeginPVTime.begin();
0621 
0622     if (beginLumiOfPVFit_ == 0)
0623       beginLumiOfPVFit_ = bpv->second;  //new begin time after removing the LS
0624     if (beginLumiOfBSFit_ == 0)
0625       beginLumiOfBSFit_ = bbs->second;
0626     if (refBStime[0] == 0)
0627       refBStime[0] = bbst->second;
0628     if (refPVtime[0] == 0)
0629       refPVtime[0] = bpvt->second;
0630 
0631   }  //same logic for average fit as above commented line
0632 
0633   map<int, std::time_t>::iterator nbbst = mapBeginBSTime.begin();
0634   map<int, std::time_t>::iterator nbpvt = mapBeginPVTime.begin();
0635 
0636   if (onlineMode_ && (nthlumi < nextlumi_))
0637     return;
0638 
0639   if (onlineMode_) {
0640     if (nthlumi > nextlumi_) {
0641       if (countLumi_ != 0 && processed_)
0642         FitAndFill(lumiSeg, lastlumi_, nextlumi_, nthlumi);
0643       nextlumi_ = nthlumi;
0644       edm::LogInfo("BeamMonitor") << "beginLuminosityBlock:: Next Lumi to Fit: " << nextlumi_ << endl;
0645       if ((StartAverage_) && refBStime[0] == 0)
0646         refBStime[0] = nbbst->second;
0647       if ((StartAverage_) && refPVtime[0] == 0)
0648         refPVtime[0] = nbpvt->second;
0649     }
0650   } else {
0651     if (processed_)
0652       FitAndFill(lumiSeg, lastlumi_, nextlumi_, nthlumi);
0653     nextlumi_ = nthlumi;
0654     edm::LogInfo("BeamMonitor") << " beginLuminosityBlock:: Next Lumi to Fit: " << nextlumi_ << endl;
0655     if ((StartAverage_) && refBStime[0] == 0)
0656       refBStime[0] = nbbst->second;
0657     if ((StartAverage_) && refPVtime[0] == 0)
0658       refPVtime[0] = nbpvt->second;
0659   }
0660 
0661   //countLumi_++;
0662   if (processed_)
0663     processed_ = false;
0664   edm::LogInfo("BeamMonitor") << " beginLuminosityBlock::  Begin of Lumi: " << nthlumi << endl;
0665 }
0666 
0667 // ----------------------------------------------------------
0668 void BeamMonitor::analyze(const Event& iEvent, const EventSetup& iSetup) {
0669   const TCDSRecord& tcdsData = iEvent.get(tcdsToken_);
0670   int beamMode = tcdsData.getBST().getBeamMode();
0671 
0672   // Check that the beamMode information is available in the event content
0673   if (beamMode == BSTRecord::BeamMode::NOMODE)
0674     edm::LogError("BeamMonitor") << "No BeamMode identified from BSTRecord!"
0675                                     "Please check that the event content has the raw data from TCDS FEDs (1024,1025)!";
0676 
0677   // Check if stable beams are present
0678   if (beamMode == BSTRecord::BeamMode::STABLE)
0679     logToDb_ = true;
0680 
0681   const int nthlumi = iEvent.luminosityBlock();
0682   if (onlineMode_ && (nthlumi < nextlumi_)) {
0683     edm::LogInfo("BeamMonitor") << "analyze::  Spilt event from previous lumi section!" << std::endl;
0684     return;
0685   }
0686   if (onlineMode_ && (nthlumi > nextlumi_)) {
0687     edm::LogInfo("BeamMonitor") << "analyze::  Spilt event from next lumi section!!!" << std::endl;
0688     return;
0689   }
0690 
0691   countEvt_++;
0692   theBeamFitter->readEvent(
0693       iEvent);  //Remember when track fitter read the event in the same place the PVFitter read the events !!!!!!!!!
0694 
0695   Handle<reco::BeamSpot> recoBeamSpotHandle;
0696   iEvent.getByToken(bsSrc_, recoBeamSpotHandle);
0697   refBS = *recoBeamSpotHandle;
0698 
0699   //------Cut Flow Table filled every event!--------------------------------------
0700   {
0701     // Make a copy of the cut flow table from the beam fitter.
0702     auto tmphisto = static_cast<TH1F*>(theBeamFitter->getCutFlow());
0703     cutFlowTable->getTH1()->SetBins(
0704         tmphisto->GetNbinsX(), tmphisto->GetXaxis()->GetXmin(), tmphisto->GetXaxis()->GetXmax());
0705     // Update the bin labels
0706     if (countEvt_ == 1)  // SetLabel just once
0707       for (int n = 0; n < tmphisto->GetNbinsX(); n++)
0708         cutFlowTable->setBinLabel(n + 1, tmphisto->GetXaxis()->GetBinLabel(n + 1), 1);
0709     cutFlowTable->Reset();
0710     cutFlowTable->getTH1()->Add(tmphisto);
0711   }
0712 
0713   //----Reco tracks -------------------------------------
0714   Handle<reco::TrackCollection> TrackCollection;
0715   iEvent.getByToken(tracksLabel_, TrackCollection);
0716   const reco::TrackCollection* tracks = TrackCollection.product();
0717   for (reco::TrackCollection::const_iterator track = tracks->begin(); track != tracks->end(); ++track) {
0718     h_trkPt->Fill(track->pt());  //no need to change  here for average bs
0719     h_trkVz->Fill(track->vz());
0720   }
0721 
0722   //-------HLT Trigger --------------------------------
0723   edm::Handle<TriggerResults> triggerResults;
0724   bool JetTrigPass = false;
0725   if (iEvent.getByToken(hltSrc_, triggerResults)) {
0726     const edm::TriggerNames& trigNames = iEvent.triggerNames(*triggerResults);
0727     for (unsigned int i = 0; i < triggerResults->size(); i++) {
0728       const std::string& trigName = trigNames.triggerName(i);
0729 
0730       if (JetTrigPass)
0731         continue;
0732 
0733       for (size_t t = 0; t < jetTrigger_.size(); ++t) {
0734         if (JetTrigPass)
0735           continue;
0736 
0737         string string_search(jetTrigger_[t]);
0738         size_t found = trigName.find(string_search);
0739 
0740         if (found != string::npos) {
0741           int thisTrigger_ = trigNames.triggerIndex(trigName);
0742           if (triggerResults->accept(thisTrigger_))
0743             JetTrigPass = true;
0744         }  //if trigger found
0745       }  //for(t=0;..)
0746     }  //for(i=0; ..)
0747   }  //if trigger colleciton exist)
0748 
0749   //------ Primary Vertices-------
0750   edm::Handle<reco::VertexCollection> PVCollection;
0751 
0752   if (iEvent.getByToken(pvSrc_, PVCollection)) {
0753     int nPVcount = 0;
0754     int nPVcount_ST = 0;  //For Single Trigger(hence ST)
0755 
0756     for (reco::VertexCollection::const_iterator pv = PVCollection->begin(); pv != PVCollection->end(); ++pv) {
0757       //--- vertex selection
0758       if (pv->isFake() || pv->tracksSize() == 0)
0759         continue;
0760       nPVcount++;  // count non fake pv:
0761 
0762       if (JetTrigPass)
0763         nPVcount_ST++;  //non-fake pv with a specific trigger
0764 
0765       if (pv->ndof() < minVtxNdf_ || (pv->ndof() + 3.) / pv->tracksSize() < 2 * minVtxWgt_)
0766         continue;
0767 
0768       //Fill this map to store xyx for pv so that later we can remove the first one for run aver
0769       mapPVx[countLumi_].push_back(pv->x());
0770       mapPVy[countLumi_].push_back(pv->y());
0771       mapPVz[countLumi_].push_back(pv->z());
0772 
0773       if (!StartAverage_) {  //for first N LS
0774         h_PVx[0]->Fill(pv->x());
0775         h_PVy[0]->Fill(pv->y());
0776         h_PVz[0]->Fill(pv->z());
0777         h_PVxz->Fill(pv->z(), pv->x());
0778         h_PVyz->Fill(pv->z(), pv->y());
0779       }  //for first N LiS
0780       else {
0781         h_PVxz->Fill(pv->z(), pv->x());
0782         h_PVyz->Fill(pv->z(), pv->y());
0783       }
0784 
0785     }  //loop over pvs
0786 
0787     h_nVtx->Fill(nPVcount * 1.);  //no need to change it for average BS
0788 
0789     mapNPV[countLumi_].push_back((nPVcount_ST));
0790 
0791     if (!StartAverage_) {
0792       h_nVtx_st->Fill(nPVcount_ST * 1.);
0793     }
0794 
0795   }  //if pv collection is availaable
0796 
0797   if (StartAverage_) {
0798     map<int, std::vector<float> >::iterator itpvx = mapPVx.begin();
0799     map<int, std::vector<float> >::iterator itpvy = mapPVy.begin();
0800     map<int, std::vector<float> >::iterator itpvz = mapPVz.begin();
0801 
0802     map<int, std::vector<int> >::iterator itbspvinfo = mapNPV.begin();
0803 
0804     if ((int)mapPVx.size() > resetFitNLumi_) {  //sometimes the events is not there but LS is there!
0805       mapPVx.erase(itpvx);
0806       mapPVy.erase(itpvy);
0807       mapPVz.erase(itpvz);
0808       mapNPV.erase(itbspvinfo);
0809     }  //loop over Last N lumi collected
0810 
0811   }  //StartAverage==true
0812 
0813   processed_ = true;
0814 }
0815 
0816 //--------------------------------------------------------
0817 void BeamMonitor::endLuminosityBlock(const LuminosityBlock& lumiSeg, const EventSetup& iSetup) {
0818   int nthlumi = lumiSeg.id().luminosityBlock();
0819   edm::LogInfo("BeamMonitor") << "endLuminosityBlock:: Lumi of the last event before endLuminosityBlock: " << nthlumi
0820                               << endl;
0821 
0822   if (onlineMode_ && nthlumi < nextlumi_)
0823     return;
0824   const edm::TimeValue_t fendtimestamp = lumiSeg.endTime().value();
0825   const std::time_t fendtime = fendtimestamp >> 32;
0826   tmpTime = refBStime[1] = refPVtime[1] = fendtime;
0827 
0828   // end DB logger
0829   if (onlineDbService_.isAvailable() && logToDb_ && loggerActive_) {
0830     onlineDbService_->logger().logInfo() << "BeamMonitor::endLuminosityBlock";
0831     onlineDbService_->logger().end(DBloggerReturn_);
0832   }
0833 }
0834 
0835 //--------------------------------------------------------
0836 void BeamMonitor::FitAndFill(const LuminosityBlock& lumiSeg, int& lastlumi, int& nextlumi, int& nthlumi) {
0837   if (onlineMode_ && (nthlumi <= nextlumi))
0838     return;
0839 
0840   //set the correct run number when no event in the LS for fake output
0841   if ((processed_) && theBeamFitter->getRunNumber() != frun)
0842     theBeamFitter->setRun(frun);
0843 
0844   int currentlumi = nextlumi;
0845   edm::LogInfo("BeamMonitor") << "FitAndFill::  Lumi of the current fit: " << currentlumi << endl;
0846   lastlumi = currentlumi;
0847   endLumiOfBSFit_ = currentlumi;
0848   endLumiOfPVFit_ = currentlumi;
0849 
0850   //---------Fix for Runninv average-------------
0851   mapLSPVStoreSize[countLumi_] = theBeamFitter->getPVvectorSize();
0852 
0853   if (StartAverage_) {
0854     std::map<int, std::size_t>::iterator rmLSPVi = mapLSPVStoreSize.begin();
0855     size_t SizeToRemovePV = rmLSPVi->second;
0856     for (std::map<int, std::size_t>::iterator rmLSPVe = mapLSPVStoreSize.end(); ++rmLSPVi != rmLSPVe;)
0857       rmLSPVi->second -= SizeToRemovePV;
0858 
0859     theBeamFitter->resizePVvector(SizeToRemovePV);
0860 
0861     map<int, std::size_t>::iterator tmpItpv = mapLSPVStoreSize.begin();
0862     mapLSPVStoreSize.erase(tmpItpv);
0863   }
0864   if (debug_)
0865     edm::LogInfo("BeamMonitor") << "FitAndFill:: Size of thePVvector After removing the PVs = "
0866                                 << theBeamFitter->getPVvectorSize() << endl;
0867 
0868   //lets filt the PV for GUI here: It was in analyzer in preivous versiton but moved here due to absence of event in some lumis, works OK
0869   bool resetHistoFlag_ = false;
0870   if ((int)mapPVx.size() >= resetFitNLumi_ && (StartAverage_)) {
0871     h_PVx[0]->Reset();
0872     h_PVy[0]->Reset();
0873     h_PVz[0]->Reset();
0874     h_nVtx_st->Reset();
0875     resetHistoFlag_ = true;
0876   }
0877 
0878   int MaxPVs = 0;
0879   int countEvtLastNLS_ = 0;
0880   int countTotPV_ = 0;
0881 
0882   std::map<int, std::vector<int> >::iterator mnpv = mapNPV.begin();
0883   std::map<int, std::vector<float> >::iterator mpv2 = mapPVy.begin();
0884   std::map<int, std::vector<float> >::iterator mpv3 = mapPVz.begin();
0885 
0886   for (std::map<int, std::vector<float> >::iterator mpv1 = mapPVx.begin(); mpv1 != mapPVx.end();
0887        ++mpv1, ++mpv2, ++mpv3, ++mnpv) {
0888     std::vector<float>::iterator mpvs2 = (mpv2->second).begin();
0889     std::vector<float>::iterator mpvs3 = (mpv3->second).begin();
0890     for (std::vector<float>::iterator mpvs1 = (mpv1->second).begin(); mpvs1 != (mpv1->second).end();
0891          ++mpvs1, ++mpvs2, ++mpvs3) {
0892       if (resetHistoFlag_) {
0893         h_PVx[0]->Fill(*mpvs1);  //these histogram are reset after StartAverage_ flag is ON
0894         h_PVy[0]->Fill(*mpvs2);
0895         h_PVz[0]->Fill(*mpvs3);
0896       }
0897     }  //loop over second
0898 
0899     //Do the same here for nPV distr.
0900     for (std::vector<int>::iterator mnpvs = (mnpv->second).begin(); mnpvs != (mnpv->second).end(); ++mnpvs) {
0901       if ((*mnpvs > 0) && (resetHistoFlag_))
0902         h_nVtx_st->Fill((*mnpvs) * (1.0));
0903       countEvtLastNLS_++;
0904       countTotPV_ += (*mnpvs);
0905       if ((*mnpvs) > MaxPVs)
0906         MaxPVs = (*mnpvs);
0907     }  //loop over second of mapNPV
0908 
0909   }  //loop over last N lumis
0910 
0911   char tmpTitlePV[100];
0912   sprintf(tmpTitlePV, "%s %i %s %i", "Num. of reco. vertices for LS: ", beginLumiOfPVFit_, " to ", endLumiOfPVFit_);
0913   h_nVtx_st->setAxisTitle(tmpTitlePV, 1);
0914 
0915   std::vector<float> DipPVInfo_;
0916   DipPVInfo_.clear();
0917 
0918   if (countTotPV_ != 0) {
0919     DipPVInfo_.push_back((float)countEvtLastNLS_);
0920     DipPVInfo_.push_back(h_nVtx_st->getMean());
0921     DipPVInfo_.push_back(h_nVtx_st->getMeanError());
0922     DipPVInfo_.push_back(h_nVtx_st->getRMS());
0923     DipPVInfo_.push_back(h_nVtx_st->getRMSError());
0924     DipPVInfo_.push_back((float)MaxPVs);
0925     DipPVInfo_.push_back((float)countTotPV_);
0926   } else {
0927     for (size_t i = 0; i < 7; i++) {
0928       if (i > 0) {
0929         DipPVInfo_.push_back(0.);
0930       } else {
0931         DipPVInfo_.push_back((float)countEvtLastNLS_);
0932       }
0933     }
0934   }
0935   theBeamFitter->SetPVInfo(DipPVInfo_);
0936 
0937   if (onlineMode_) {  // filling LS gap
0938     // FIXME: need to add protection for the case if the gap is at the resetting LS!
0939     const int countLS_bs = hs[k_x0_lumi]->getTH1()->GetEntries();
0940     const int countLS_pv = hs[k_PVx_lumi]->getTH1()->GetEntries();
0941     edm::LogInfo("BeamMonitor") << "FitAndFill:: countLS_bs = " << countLS_bs << " ; countLS_pv = " << countLS_pv
0942                                 << std::endl;
0943     int LSgap_bs = currentlumi / fitNLumi_ - countLS_bs;
0944     int LSgap_pv = currentlumi / fitPVNLumi_ - countLS_pv;
0945     if (currentlumi % fitNLumi_ == 0)
0946       LSgap_bs--;
0947     if (currentlumi % fitPVNLumi_ == 0)
0948       LSgap_pv--;
0949     edm::LogInfo("BeamMonitor") << "FitAndFill::  LSgap_bs = " << LSgap_bs << " ; LSgap_pv = " << LSgap_pv << std::endl;
0950     // filling previous fits if LS gap ever exists
0951     for (int ig = 0; ig < LSgap_bs; ig++) {
0952       hs[k_x0_lumi]->ShiftFillLast(0., 0., fitNLumi_);  //x0 , x0err, fitNLumi_;  see DQMCore....
0953       hs[k_y0_lumi]->ShiftFillLast(0., 0., fitNLumi_);
0954       hs[k_z0_lumi]->ShiftFillLast(0., 0., fitNLumi_);
0955       hs[k_sigmaX0_lumi]->ShiftFillLast(0., 0., fitNLumi_);
0956       hs[k_sigmaY0_lumi]->ShiftFillLast(0., 0., fitNLumi_);
0957       hs[k_sigmaZ0_lumi]->ShiftFillLast(0., 0., fitNLumi_);
0958       h_nVtx_lumi->ShiftFillLast(0., 0., fitNLumi_);
0959     }
0960     for (int ig = 0; ig < LSgap_pv; ig++) {
0961       hs[k_PVx_lumi]->ShiftFillLast(0., 0., fitPVNLumi_);
0962       hs[k_PVy_lumi]->ShiftFillLast(0., 0., fitPVNLumi_);
0963       hs[k_PVz_lumi]->ShiftFillLast(0., 0., fitPVNLumi_);
0964     }
0965     const int previousLS = h_nTrk_lumi->getTH1()->GetEntries();
0966     for (int i = 1; i < (currentlumi - previousLS);
0967          i++)  //if (current-previoius)= 1 then never go inside the for loop!!!!!!!!!!!
0968       h_nTrk_lumi->ShiftFillLast(nthBSTrk_);
0969   }
0970 
0971   edm::LogInfo("BeamMonitor") << "FitAndFill:: Time lapsed since last scroll = " << tmpTime - refTime << std::endl;
0972 
0973   if (testScroll(tmpTime, refTime)) {
0974     scrollTH1(hs[k_x0_time]->getTH1(), refTime);
0975     scrollTH1(hs[k_y0_time]->getTH1(), refTime);
0976     scrollTH1(hs[k_z0_time]->getTH1(), refTime);
0977     scrollTH1(hs[k_sigmaX0_time]->getTH1(), refTime);
0978     scrollTH1(hs[k_sigmaY0_time]->getTH1(), refTime);
0979     scrollTH1(hs[k_sigmaZ0_time]->getTH1(), refTime);
0980     scrollTH1(hs[k_PVx_time]->getTH1(), refTime);
0981     scrollTH1(hs[k_PVy_time]->getTH1(), refTime);
0982     scrollTH1(hs[k_PVz_time]->getTH1(), refTime);
0983   }
0984 
0985   bool doPVFit = false;
0986 
0987   if (fitPVNLumi_ > 0) {
0988     if (onlineMode_) {
0989       if (currentlumi % fitPVNLumi_ == 0)
0990         doPVFit = true;
0991     } else if (countLumi_ % fitPVNLumi_ == 0)
0992       doPVFit = true;
0993   } else
0994     doPVFit = true;
0995 
0996   if (doPVFit) {
0997     edm::LogInfo("BeamMonitor") << "FitAndFill:: Do PV Fitting for LS = " << beginLumiOfPVFit_ << " to "
0998                                 << endLumiOfPVFit_ << std::endl;
0999     // Primary Vertex Fit:
1000     if (h_PVx[0]->getTH1()->GetEntries() > minNrVertices_) {
1001       pvResults->Reset();
1002       char tmpTitle[50];
1003       sprintf(
1004           tmpTitle, "%s %i %s %i", "Fitted Primary Vertex (cm) of LS: ", beginLumiOfPVFit_, " to ", endLumiOfPVFit_);
1005       pvResults->setAxisTitle(tmpTitle, 1);
1006 
1007       std::unique_ptr<TF1> fgaus{new TF1("fgaus", "gaus")};
1008       double mean, width, meanErr, widthErr;
1009       fgaus->SetLineColor(4);
1010       h_PVx[0]->getTH1()->Fit(fgaus.get(), "QLM0");
1011       mean = fgaus->GetParameter(1);
1012       width = fgaus->GetParameter(2);
1013       meanErr = fgaus->GetParError(1);
1014       widthErr = fgaus->GetParError(2);
1015 
1016       hs[k_PVx_lumi]->ShiftFillLast(mean, width, fitPVNLumi_);
1017       hs[k_PVx_lumi_all]->setBinContent(currentlumi, mean);
1018       hs[k_PVx_lumi_all]->setBinError(currentlumi, width);
1019       int nthBin = tmpTime - refTime;
1020       if (nthBin < 0)
1021         edm::LogInfo("BeamMonitor") << "FitAndFill::  Event time outside current range of time histograms!"
1022                                     << std::endl;
1023       if (nthBin > 0) {
1024         hs[k_PVx_time]->setBinContent(nthBin, mean);
1025         hs[k_PVx_time]->setBinError(nthBin, width);
1026       }
1027       int jthBin = tmpTime - startTime;
1028       if (jthBin > 0) {
1029         hs[k_PVx_time_all]->setBinContent(jthBin, mean);
1030         hs[k_PVx_time_all]->setBinError(jthBin, width);
1031       }
1032       pvResults->setBinContent(1, 6, mean);
1033       pvResults->setBinContent(1, 3, width);
1034       pvResults->setBinContent(2, 6, meanErr);
1035       pvResults->setBinContent(2, 3, widthErr);
1036 
1037       {
1038         // snap shot of the fit
1039         auto tmphisto = h_PVx[0]->getTH1F();
1040         h_PVx[1]->getTH1()->SetBins(
1041             tmphisto->GetNbinsX(), tmphisto->GetXaxis()->GetXmin(), tmphisto->GetXaxis()->GetXmax());
1042         h_PVx[1]->Reset();
1043         h_PVx[1]->getTH1()->Add(tmphisto);
1044         h_PVx[1]->getTH1()->Fit(fgaus.get(), "QLM");
1045       }
1046 
1047       h_PVy[0]->getTH1()->Fit(fgaus.get(), "QLM0");
1048       mean = fgaus->GetParameter(1);
1049       width = fgaus->GetParameter(2);
1050       meanErr = fgaus->GetParError(1);
1051       widthErr = fgaus->GetParError(2);
1052       hs[k_PVy_lumi]->ShiftFillLast(mean, width, fitPVNLumi_);
1053       hs[k_PVy_lumi_all]->setBinContent(currentlumi, mean);
1054       hs[k_PVy_lumi_all]->setBinError(currentlumi, width);
1055       if (nthBin > 0) {
1056         hs[k_PVy_time]->setBinContent(nthBin, mean);
1057         hs[k_PVy_time]->setBinError(nthBin, width);
1058       }
1059       if (jthBin > 0) {
1060         hs[k_PVy_time_all]->setBinContent(jthBin, mean);
1061         hs[k_PVy_time_all]->setBinError(jthBin, width);
1062       }
1063       pvResults->setBinContent(1, 5, mean);
1064       pvResults->setBinContent(1, 2, width);
1065       pvResults->setBinContent(2, 5, meanErr);
1066       pvResults->setBinContent(2, 2, widthErr);
1067       // snap shot of the fit
1068       {
1069         auto tmphisto = h_PVy[0]->getTH1F();
1070         h_PVy[1]->getTH1()->SetBins(
1071             tmphisto->GetNbinsX(), tmphisto->GetXaxis()->GetXmin(), tmphisto->GetXaxis()->GetXmax());
1072         h_PVy[1]->Reset();
1073         h_PVy[1]->getTH1()->Add(tmphisto);
1074         h_PVy[1]->getTH1()->Fit(fgaus.get(), "QLM");
1075       }
1076 
1077       h_PVz[0]->getTH1()->Fit(fgaus.get(), "QLM0");
1078       mean = fgaus->GetParameter(1);
1079       width = fgaus->GetParameter(2);
1080       meanErr = fgaus->GetParError(1);
1081       widthErr = fgaus->GetParError(2);
1082       hs[k_PVz_lumi]->ShiftFillLast(mean, width, fitPVNLumi_);
1083       hs[k_PVz_lumi_all]->setBinContent(currentlumi, mean);
1084       hs[k_PVz_lumi_all]->setBinError(currentlumi, width);
1085       if (nthBin > 0) {
1086         hs[k_PVz_time]->setBinContent(nthBin, mean);
1087         hs[k_PVz_time]->setBinError(nthBin, width);
1088       }
1089       if (jthBin > 0) {
1090         hs[k_PVz_time_all]->setBinContent(jthBin, mean);
1091         hs[k_PVz_time_all]->setBinError(jthBin, width);
1092       }
1093       pvResults->setBinContent(1, 4, mean);
1094       pvResults->setBinContent(1, 1, width);
1095       pvResults->setBinContent(2, 4, meanErr);
1096       pvResults->setBinContent(2, 1, widthErr);
1097       {
1098         // snap shot of the fit
1099         auto tmphisto = h_PVz[0]->getTH1F();
1100         h_PVz[1]->getTH1()->SetBins(
1101             tmphisto->GetNbinsX(), tmphisto->GetXaxis()->GetXmin(), tmphisto->GetXaxis()->GetXmax());
1102         h_PVz[1]->Reset();
1103         h_PVz[1]->getTH1()->Add(tmphisto);
1104         h_PVz[1]->getTH1()->Fit(fgaus.get(), "QLM");
1105       }
1106     }  //check if found min Vertices
1107   }  //do PVfit
1108 
1109   if ((resetPVNLumi_ > 0 && countLumi_ == resetPVNLumi_) || StartAverage_) {
1110     beginLumiOfPVFit_ = 0;
1111     refPVtime[0] = 0;
1112   }
1113 
1114   //---------Readjustment of theBSvector, RefTime, beginLSofFit---------
1115   vector<BSTrkParameters> theBSvector1 = theBeamFitter->getBSvector();
1116   mapLSBSTrkSize[countLumi_] = (theBSvector1.size());
1117   size_t PreviousRecords = 0;  //needed to fill nth record of tracks in GUI
1118 
1119   if (StartAverage_) {
1120     size_t SizeToRemove = 0;
1121     std::map<int, std::size_t>::iterator rmls = mapLSBSTrkSize.begin();
1122     SizeToRemove = rmls->second;
1123     if (debug_)
1124       edm::LogInfo("BeamMonitor") << "  The size to remove is =  " << SizeToRemove << endl;
1125     int changedAfterThis = 0;
1126     for (std::map<int, std::size_t>::iterator rmLS = mapLSBSTrkSize.begin(); rmLS != mapLSBSTrkSize.end();
1127          ++rmLS, ++changedAfterThis) {
1128       if (changedAfterThis > 0) {
1129         (rmLS->second) = (rmLS->second) - SizeToRemove;
1130         if ((mapLSBSTrkSize.size() - (size_t)changedAfterThis) == 2)
1131           PreviousRecords = (rmLS->second);
1132       }
1133     }
1134 
1135     theBeamFitter->resizeBSvector(SizeToRemove);
1136 
1137     map<int, std::size_t>::iterator tmpIt = mapLSBSTrkSize.begin();
1138     mapLSBSTrkSize.erase(tmpIt);
1139 
1140     std::pair<int, int> checkfitLS = theBeamFitter->getFitLSRange();
1141     std::pair<time_t, time_t> checkfitTime = theBeamFitter->getRefTime();
1142     theBeamFitter->setFitLSRange(beginLumiOfBSFit_, checkfitLS.second);
1143     theBeamFitter->setRefTime(refBStime[0], checkfitTime.second);
1144   }
1145 
1146   //Fill the track for this fit
1147   vector<BSTrkParameters> theBSvector = theBeamFitter->getBSvector();
1148   h_nTrk_lumi->ShiftFillLast(theBSvector.size());
1149 
1150   if (debug_)
1151     edm::LogInfo("BeamMonitor") << "FitAndFill::   Size of  theBSViector.size()  After =" << theBSvector.size() << endl;
1152 
1153   bool countFitting = false;
1154   if (theBSvector.size() >= PreviousRecords && theBSvector.size() >= min_Ntrks_) {
1155     countFitting = true;
1156   }
1157 
1158   //---Fix for Cut Flow Table for Running average in a same way//the previous code  has problem for resetting!!!
1159   mapLSCF[countLumi_] = *theBeamFitter->getCutFlow();
1160   if (StartAverage_ && !mapLSCF.empty()) {
1161     const TH1F& cutFlowToSubtract = mapLSCF.begin()->second;
1162     // Subtract the last cut flow from all of the others.
1163     std::map<int, TH1F>::iterator cf = mapLSCF.begin();
1164     // Start on second entry
1165     for (; cf != mapLSCF.end(); ++cf) {
1166       cf->second.Add(&cutFlowToSubtract, -1);
1167     }
1168     theBeamFitter->subtractFromCutFlow(&cutFlowToSubtract);
1169     // Remove the obsolete lumi section
1170     mapLSCF.erase(mapLSCF.begin());
1171   }
1172 
1173   if (resetHistos_) {
1174     h_d0_phi0->Reset();
1175     h_vx_vy->Reset();
1176     h_vx_dz->Reset();
1177     h_vy_dz->Reset();
1178     h_trk_z0->Reset();
1179     resetHistos_ = false;
1180   }
1181 
1182   if (StartAverage_)
1183     nthBSTrk_ = PreviousRecords;  //after average proccess is ON//for 2-6 LS fit PreviousRecords is size from 2-5 LS
1184 
1185   edm::LogInfo("BeamMonitor") << " The Previous Recored for this fit is  =" << nthBSTrk_ << endl;
1186 
1187   unsigned int itrk = 0;
1188   for (vector<BSTrkParameters>::const_iterator BSTrk = theBSvector.begin(); BSTrk != theBSvector.end();
1189        ++BSTrk, ++itrk) {
1190     if (itrk >= nthBSTrk_) {  //fill for this record only !!
1191       h_d0_phi0->Fill(BSTrk->phi0(), BSTrk->d0());
1192       double vx = BSTrk->vx();
1193       double vy = BSTrk->vy();
1194       double z0 = BSTrk->z0();
1195       h_vx_vy->Fill(vx, vy);
1196       h_vx_dz->Fill(z0, vx);
1197       h_vy_dz->Fill(z0, vy);
1198       h_trk_z0->Fill(z0);
1199     }
1200   }
1201 
1202   nthBSTrk_ = theBSvector.size();  // keep track of num of tracks filled so far
1203 
1204   edm::LogInfo("BeamMonitor") << " The Current Recored for this fit is  =" << nthBSTrk_ << endl;
1205 
1206   if (countFitting)
1207     edm::LogInfo("BeamMonitor") << "FitAndFill::  Num of tracks collected = " << nthBSTrk_ << endl;
1208 
1209   if (fitNLumi_ > 0) {
1210     if (onlineMode_) {
1211       if (currentlumi % fitNLumi_ != 0) {
1212         //  for (std::map<TString,MonitorElement*>::iterator itAll = hs.begin();
1213         //       itAll != hs.end(); ++itAll) {
1214         //    if ((*itAll).first.Contains("all")) {
1215         //      (*itAll).second->setBinContent(currentlumi,0.);
1216         //      (*itAll).second->setBinError(currentlumi,0.);
1217         //    }
1218         //  }
1219         return;
1220       }
1221     } else if (countLumi_ % fitNLumi_ != 0)
1222       return;
1223   }
1224 
1225   edm::LogInfo("BeamMonitor") << "FitAndFill:: [DebugTime] refBStime[0] = " << refBStime[0]
1226                               << "; address =  " << &refBStime[0] << std::endl;
1227   edm::LogInfo("BeamMonitor") << "FitAndFill:: [DebugTime] refBStime[1] = " << refBStime[1]
1228                               << "; address =  " << &refBStime[1] << std::endl;
1229 
1230   //Fill for all LS even if fit fails
1231   h_nVtx_lumi->ShiftFillLast((theBeamFitter->getPVvectorSize()), 0., fitNLumi_);
1232   h_nVtx_lumi_all->setBinContent(currentlumi, (theBeamFitter->getPVvectorSize()));
1233 
1234   if (countFitting) {
1235     nFits_++;
1236     std::pair<int, int> fitLS = theBeamFitter->getFitLSRange();
1237     edm::LogInfo("BeamMonitor") << "FitAndFill::  [BeamFitter] Do BeamSpot Fit for LS = " << fitLS.first << " to "
1238                                 << fitLS.second << std::endl;
1239     edm::LogInfo("BeamMonitor") << "FitAndFill::  [BeamMonitor] Do BeamSpot Fit for LS = " << beginLumiOfBSFit_
1240                                 << " to " << endLumiOfBSFit_ << std::endl;
1241 
1242     //Now Run the PV and Track Fitter over the collected tracks and pvs
1243     if (theBeamFitter->runPVandTrkFitter()) {
1244       reco::BeamSpot bs = theBeamFitter->getBeamSpot();
1245       if (bs.type() > 0)  // with good beamwidth fit
1246         preBS = bs;       // cache good fit results
1247 
1248       edm::LogInfo("BeamMonitor") << "\n RESULTS OF DEFAULT FIT:" << endl;
1249       edm::LogInfo("BeamMonitor") << bs << endl;
1250       edm::LogInfo("BeamMonitor") << "[BeamFitter] fitting done \n" << endl;
1251 
1252       hs[k_x0_lumi]->ShiftFillLast(bs.x0(), bs.x0Error(), fitNLumi_);
1253       hs[k_y0_lumi]->ShiftFillLast(bs.y0(), bs.y0Error(), fitNLumi_);
1254       hs[k_z0_lumi]->ShiftFillLast(bs.z0(), bs.z0Error(), fitNLumi_);
1255       hs[k_sigmaX0_lumi]->ShiftFillLast(bs.BeamWidthX(), bs.BeamWidthXError(), fitNLumi_);
1256       hs[k_sigmaY0_lumi]->ShiftFillLast(bs.BeamWidthY(), bs.BeamWidthYError(), fitNLumi_);
1257       hs[k_sigmaZ0_lumi]->ShiftFillLast(bs.sigmaZ(), bs.sigmaZ0Error(), fitNLumi_);
1258       hs[k_x0_lumi_all]->setBinContent(currentlumi, bs.x0());
1259       hs[k_x0_lumi_all]->setBinError(currentlumi, bs.x0Error());
1260       hs[k_y0_lumi_all]->setBinContent(currentlumi, bs.y0());
1261       hs[k_y0_lumi_all]->setBinError(currentlumi, bs.y0Error());
1262       hs[k_z0_lumi_all]->setBinContent(currentlumi, bs.z0());
1263       hs[k_z0_lumi_all]->setBinError(currentlumi, bs.z0Error());
1264       hs[k_sigmaX0_lumi_all]->setBinContent(currentlumi, bs.BeamWidthX());
1265       hs[k_sigmaX0_lumi_all]->setBinError(currentlumi, bs.BeamWidthXError());
1266       hs[k_sigmaY0_lumi_all]->setBinContent(currentlumi, bs.BeamWidthY());
1267       hs[k_sigmaY0_lumi_all]->setBinError(currentlumi, bs.BeamWidthYError());
1268       hs[k_sigmaZ0_lumi_all]->setBinContent(currentlumi, bs.sigmaZ());
1269       hs[k_sigmaZ0_lumi_all]->setBinError(currentlumi, bs.sigmaZ0Error());
1270 
1271       int nthBin = tmpTime - refTime;
1272       if (nthBin > 0) {
1273         hs[k_x0_time]->setBinContent(nthBin, bs.x0());
1274         hs[k_y0_time]->setBinContent(nthBin, bs.y0());
1275         hs[k_z0_time]->setBinContent(nthBin, bs.z0());
1276         hs[k_sigmaX0_time]->setBinContent(nthBin, bs.BeamWidthX());
1277         hs[k_sigmaY0_time]->setBinContent(nthBin, bs.BeamWidthY());
1278         hs[k_sigmaZ0_time]->setBinContent(nthBin, bs.sigmaZ());
1279         hs[k_x0_time]->setBinError(nthBin, bs.x0Error());
1280         hs[k_y0_time]->setBinError(nthBin, bs.y0Error());
1281         hs[k_z0_time]->setBinError(nthBin, bs.z0Error());
1282         hs[k_sigmaX0_time]->setBinError(nthBin, bs.BeamWidthXError());
1283         hs[k_sigmaY0_time]->setBinError(nthBin, bs.BeamWidthYError());
1284         hs[k_sigmaZ0_time]->setBinError(nthBin, bs.sigmaZ0Error());
1285       }
1286 
1287       int jthBin = tmpTime - startTime;
1288       if (jthBin > 0) {
1289         hs[k_x0_time_all]->setBinContent(jthBin, bs.x0());
1290         hs[k_y0_time_all]->setBinContent(jthBin, bs.y0());
1291         hs[k_z0_time_all]->setBinContent(jthBin, bs.z0());
1292         hs[k_sigmaX0_time_all]->setBinContent(jthBin, bs.BeamWidthX());
1293         hs[k_sigmaY0_time_all]->setBinContent(jthBin, bs.BeamWidthY());
1294         hs[k_sigmaZ0_time_all]->setBinContent(jthBin, bs.sigmaZ());
1295         hs[k_x0_time_all]->setBinError(jthBin, bs.x0Error());
1296         hs[k_y0_time_all]->setBinError(jthBin, bs.y0Error());
1297         hs[k_z0_time_all]->setBinError(jthBin, bs.z0Error());
1298         hs[k_sigmaX0_time_all]->setBinError(jthBin, bs.BeamWidthXError());
1299         hs[k_sigmaY0_time_all]->setBinError(jthBin, bs.BeamWidthYError());
1300         hs[k_sigmaZ0_time_all]->setBinError(jthBin, bs.sigmaZ0Error());
1301       }
1302 
1303       h_x0->Fill(bs.x0());
1304       h_y0->Fill(bs.y0());
1305       h_z0->Fill(bs.z0());
1306       if (bs.type() > 0) {  // with good beamwidth fit
1307         h_sigmaX0->Fill(bs.BeamWidthX());
1308         h_sigmaY0->Fill(bs.BeamWidthY());
1309       }
1310       h_sigmaZ0->Fill(bs.sigmaZ());
1311 
1312       if (nthBSTrk_ >= 2 * min_Ntrks_) {
1313         double amp = std::sqrt(bs.x0() * bs.x0() + bs.y0() * bs.y0());
1314         double alpha = std::atan2(bs.y0(), bs.x0());
1315         std::unique_ptr<TF1> f1{new TF1("f1", "[0]*sin(x-[1])", -3.14, 3.14)};
1316         f1->SetParameters(amp, alpha);
1317         f1->SetParLimits(0, amp - 0.1, amp + 0.1);
1318         f1->SetParLimits(1, alpha - 0.577, alpha + 0.577);
1319         f1->SetLineColor(4);
1320         h_d0_phi0->getTProfile()->Fit(f1.get(), "QR");
1321 
1322         double mean = bs.z0();
1323         double width = bs.sigmaZ();
1324         std::unique_ptr<TF1> fgaus{new TF1("fgaus", "gaus")};
1325         fgaus->SetParameters(mean, width);
1326         fgaus->SetLineColor(4);
1327         h_trk_z0->getTH1()->Fit(fgaus.get(), "QLRM", "", mean - 3 * width, mean + 3 * width);
1328       }
1329 
1330       fitResults->Reset();
1331       std::pair<int, int> LSRange = theBeamFitter->getFitLSRange();
1332       char tmpTitle[50];
1333       sprintf(tmpTitle, "%s %i %s %i", "Fitted Beam Spot (cm) of LS: ", LSRange.first, " to ", LSRange.second);
1334       fitResults->setAxisTitle(tmpTitle, 1);
1335       fitResults->setBinContent(1, 8, bs.x0());
1336       fitResults->setBinContent(1, 7, bs.y0());
1337       fitResults->setBinContent(1, 6, bs.z0());
1338       fitResults->setBinContent(1, 5, bs.sigmaZ());
1339       fitResults->setBinContent(1, 4, bs.dxdz());
1340       fitResults->setBinContent(1, 3, bs.dydz());
1341       if (bs.type() > 0) {  // with good beamwidth fit
1342         fitResults->setBinContent(1, 2, bs.BeamWidthX());
1343         fitResults->setBinContent(1, 1, bs.BeamWidthY());
1344       } else {  // fill cached widths
1345         fitResults->setBinContent(1, 2, preBS.BeamWidthX());
1346         fitResults->setBinContent(1, 1, preBS.BeamWidthY());
1347       }
1348 
1349       fitResults->setBinContent(2, 8, bs.x0Error());
1350       fitResults->setBinContent(2, 7, bs.y0Error());
1351       fitResults->setBinContent(2, 6, bs.z0Error());
1352       fitResults->setBinContent(2, 5, bs.sigmaZ0Error());
1353       fitResults->setBinContent(2, 4, bs.dxdzError());
1354       fitResults->setBinContent(2, 3, bs.dydzError());
1355       if (bs.type() > 0) {  // with good beamwidth fit
1356         fitResults->setBinContent(2, 2, bs.BeamWidthXError());
1357         fitResults->setBinContent(2, 1, bs.BeamWidthYError());
1358       } else {  // fill cached width errors
1359         fitResults->setBinContent(2, 2, preBS.BeamWidthXError());
1360         fitResults->setBinContent(2, 1, preBS.BeamWidthYError());
1361       }
1362 
1363       // count good fit
1364       //     if (std::fabs(refBS.x0()-bs.x0())/bs.x0Error() < deltaSigCut_) { // disabled temporarily
1365       summaryContent_[0] += 1.;
1366       //     }
1367       //     if (std::fabs(refBS.y0()-bs.y0())/bs.y0Error() < deltaSigCut_) { // disabled temporarily
1368       summaryContent_[1] += 1.;
1369       //     }
1370       //     if (std::fabs(refBS.z0()-bs.z0())/bs.z0Error() < deltaSigCut_) { // disabled temporarily
1371       summaryContent_[2] += 1.;
1372       //     }
1373 
1374       // Create the BeamSpotOnlineObjects object
1375       BeamSpotOnlineObjects BSOnline;
1376       BSOnline.setLastAnalyzedLumi(LSRange.second);
1377       BSOnline.setLastAnalyzedRun(theBeamFitter->getRunNumber());
1378       BSOnline.setLastAnalyzedFill(0);  // To be updated with correct LHC Fill number
1379       BSOnline.setPosition(bs.x0(), bs.y0(), bs.z0());
1380       BSOnline.setSigmaZ(bs.sigmaZ());
1381       BSOnline.setBeamWidthX(bs.BeamWidthX());
1382       BSOnline.setBeamWidthY(bs.BeamWidthY());
1383       BSOnline.setBeamWidthXError(bs.BeamWidthXError());
1384       BSOnline.setBeamWidthYError(bs.BeamWidthYError());
1385       BSOnline.setdxdz(bs.dxdz());
1386       BSOnline.setdydz(bs.dydz());
1387       BSOnline.setType(bs.type());
1388       BSOnline.setEmittanceX(bs.emittanceX());
1389       BSOnline.setEmittanceY(bs.emittanceY());
1390       BSOnline.setBetaStar(bs.betaStar());
1391       for (int i = 0; i < 7; ++i) {
1392         for (int j = 0; j < 7; ++j) {
1393           BSOnline.setCovariance(i, j, bs.covariance(i, j));
1394         }
1395       }
1396       BSOnline.setNumTracks(theBeamFitter->getNTracks());
1397       BSOnline.setNumPVs(theBeamFitter->getNPVs());
1398       BSOnline.setUsedEvents((int)DipPVInfo_[0]);
1399       BSOnline.setMeanPV(DipPVInfo_[1]);
1400       BSOnline.setMeanErrorPV(DipPVInfo_[2]);
1401       BSOnline.setRmsPV(DipPVInfo_[3]);
1402       BSOnline.setRmsErrorPV(DipPVInfo_[4]);
1403       BSOnline.setMaxPVs((int)DipPVInfo_[5]);
1404       auto creationTime =
1405           std::chrono::duration_cast<std::chrono::microseconds>(std::chrono::system_clock::now().time_since_epoch())
1406               .count();
1407       BSOnline.setCreationTime(creationTime);
1408 
1409       std::pair<time_t, time_t> timeForDIP = theBeamFitter->getRefTime();
1410       BSOnline.setStartTimeStamp(timeForDIP.first);
1411       BSOnline.setStartTime(getGMTstring(timeForDIP.first));
1412       BSOnline.setEndTimeStamp(timeForDIP.second);
1413       BSOnline.setEndTime(getGMTstring(timeForDIP.second));
1414 
1415       std::string lumiRangeForDIP = std::to_string(LSRange.first) + " - " + std::to_string(LSRange.second);
1416       BSOnline.setLumiRange(lumiRangeForDIP);
1417 
1418       edm::LogInfo("BeamMonitor") << "FitAndFill::[PayloadCreation] BeamSpotOnline object created: \n" << std::endl;
1419       edm::LogInfo("BeamMonitor") << BSOnline << std::endl;
1420 
1421       // Create the payload for BeamSpotOnlineObjects object
1422       if (onlineDbService_.isAvailable() && (nAnalyzedLS_ < nLS_for_upload_ || nAnalyzedLS_ % nLS_for_upload_ == 0) &&
1423           logToDb_) {
1424         edm::LogInfo("BeamMonitor") << "FitAndFill::[PayloadCreation] onlineDbService available \n" << std::endl;
1425         onlineDbService_->logger().logInfo() << "BeamMonitor::FitAndFill - Lumi of the current fit: " << currentlumi;
1426         onlineDbService_->logger().logInfo()
1427             << "BeamMonitor::FitAndFill - Do PV Fitting for LS = " << beginLumiOfPVFit_ << " to " << endLumiOfPVFit_;
1428         onlineDbService_->logger().logInfo()
1429             << "BeamMonitor::FitAndFill - [BeamFitter] Do BeamSpot Fit for LS = " << LSRange.first << " to "
1430             << LSRange.second;
1431         onlineDbService_->logger().logInfo()
1432             << "BeamMonitor::FitAndFill - [BeamMonitor] Do BeamSpot Fit for LS = " << beginLumiOfBSFit_ << " to "
1433             << endLumiOfBSFit_;
1434         onlineDbService_->logger().logInfo() << "BeamMonitor - RESULTS OF DEFAULT FIT:";
1435         onlineDbService_->logger().logInfo() << "\n" << bs;
1436         onlineDbService_->logger().logInfo()
1437             << "BeamMonitor::FitAndFill - [PayloadCreation] BeamSpotOnline object created:";
1438         onlineDbService_->logger().logInfo() << "\n" << BSOnline;
1439         onlineDbService_->logger().logInfo() << "BeamMonitor - Additional parameters for DIP:";
1440         onlineDbService_->logger().logInfo() << "Events used in the fit: " << BSOnline.usedEvents();
1441         onlineDbService_->logger().logInfo() << "Mean PV               : " << BSOnline.meanPV();
1442         onlineDbService_->logger().logInfo() << "Mean PV Error         : " << BSOnline.meanErrorPV();
1443         onlineDbService_->logger().logInfo() << "Rms PV                : " << BSOnline.rmsPV();
1444         onlineDbService_->logger().logInfo() << "Rms PV Error          : " << BSOnline.rmsErrorPV();
1445         onlineDbService_->logger().logInfo() << "Max PVs               : " << BSOnline.maxPVs();
1446         onlineDbService_->logger().logInfo() << "StartTime             : " << BSOnline.startTime();
1447         onlineDbService_->logger().logInfo() << "StartTimeStamp        : " << BSOnline.startTimeStamp();
1448         onlineDbService_->logger().logInfo() << "EndTime               : " << BSOnline.endTime();
1449         onlineDbService_->logger().logInfo() << "EndTimeStamp          : " << BSOnline.endTimeStamp();
1450         onlineDbService_->logger().logInfo() << "BeamMonitor::FitAndFill - [PayloadCreation] onlineDbService available";
1451         onlineDbService_->logger().logInfo()
1452             << "BeamMonitor::FitAndFill - [PayloadCreation] SetCreationTime: " << creationTime
1453             << " [epoch in microseconds]";
1454         try {
1455           onlineDbService_->writeIOVForNextLumisection(BSOnline, recordName_);
1456           onlineDbService_->logger().logInfo()
1457               << "BeamMonitor::FitAndFill - [PayloadCreation] writeIOVForNextLumisection executed correctly";
1458         } catch (const std::exception& e) {
1459           onlineDbService_->logger().logError() << "BeamMonitor - Error writing record: " << recordName_
1460                                                 << " for Run: " << frun << " - Lumi: " << LSRange.second;
1461           onlineDbService_->logger().logError() << "Error is: " << e.what();
1462           onlineDbService_->logger().logError() << "RESULTS OF DEFAULT FIT WAS:";
1463           onlineDbService_->logger().logError() << "\n" << bs;
1464           DBloggerReturn_ = 2;
1465         }
1466       }
1467       edm::LogInfo("BeamMonitor") << "FitAndFill::[PayloadCreation] BeamSpotOnline payload created \n" << std::endl;
1468 
1469     }  //if (theBeamFitter->runPVandTrkFitter())
1470     else {  // beam fit fails
1471       reco::BeamSpot bs = theBeamFitter->getBeamSpot();
1472       edm::LogInfo("BeamMonitor") << "FitAndFill::   [BeamMonitor] Beam fit fails!!! \n" << endl;
1473       edm::LogInfo("BeamMonitor") << "FitAndFill::   [BeamMonitor] Output beam spot for DIP \n" << endl;
1474       edm::LogInfo("BeamMonitor") << bs << endl;
1475 
1476       if (onlineDbService_.isAvailable() && logToDb_) {
1477         onlineDbService_->logger().logInfo() << "BeamMonitor::FitAndFill - Beam fit fails!!!";
1478         onlineDbService_->logger().logInfo() << "BeamMonitor::FitAndFill - Output beam spot for DIP";
1479         onlineDbService_->logger().logInfo() << "\n" << bs;
1480         DBloggerReturn_ = 2;
1481       }
1482 
1483       hs[k_sigmaX0_lumi]->ShiftFillLast(bs.BeamWidthX(), bs.BeamWidthXError(), fitNLumi_);
1484       hs[k_sigmaY0_lumi]->ShiftFillLast(bs.BeamWidthY(), bs.BeamWidthYError(), fitNLumi_);
1485       hs[k_sigmaZ0_lumi]->ShiftFillLast(bs.sigmaZ(), bs.sigmaZ0Error(), fitNLumi_);
1486       hs[k_x0_lumi]->ShiftFillLast(bs.x0(), bs.x0Error(), fitNLumi_);
1487       hs[k_y0_lumi]->ShiftFillLast(bs.y0(), bs.y0Error(), fitNLumi_);
1488       hs[k_z0_lumi]->ShiftFillLast(bs.z0(), bs.z0Error(), fitNLumi_);
1489     }  // end of beam fit fails
1490 
1491   }  //-------- end of countFitting------------------------------------------
1492   else {  // no fit
1493     // Overwrite Fit LS and fit time when no event processed or no track selected
1494     theBeamFitter->setFitLSRange(beginLumiOfBSFit_, endLumiOfBSFit_);
1495     theBeamFitter->setRefTime(refBStime[0], refBStime[1]);
1496     if (theBeamFitter->runPVandTrkFitter()) {
1497     }  // Dump fake beam spot for DIP
1498     reco::BeamSpot bs = theBeamFitter->getBeamSpot();
1499     edm::LogInfo("BeamMonitor") << "FitAndFill::  [BeamMonitor] No fitting \n" << endl;
1500     edm::LogInfo("BeamMonitor") << "FitAndFill::  [BeamMonitor] Output fake beam spot for DIP \n" << endl;
1501     edm::LogInfo("BeamMonitor") << bs << endl;
1502 
1503     if (onlineDbService_.isAvailable() && logToDb_) {
1504       onlineDbService_->logger().logInfo() << "BeamMonitor::FitAndFill - No fitting";
1505       onlineDbService_->logger().logInfo() << "BeamMonitor::FitAndFill - Output fake beam spot for DIP";
1506       onlineDbService_->logger().logInfo() << "\n" << bs;
1507       DBloggerReturn_ = 2;
1508     }
1509 
1510     hs[k_sigmaX0_lumi]->ShiftFillLast(bs.BeamWidthX(), bs.BeamWidthXError(), fitNLumi_);
1511     hs[k_sigmaY0_lumi]->ShiftFillLast(bs.BeamWidthY(), bs.BeamWidthYError(), fitNLumi_);
1512     hs[k_sigmaZ0_lumi]->ShiftFillLast(bs.sigmaZ(), bs.sigmaZ0Error(), fitNLumi_);
1513     hs[k_x0_lumi]->ShiftFillLast(bs.x0(), bs.x0Error(), fitNLumi_);
1514     hs[k_y0_lumi]->ShiftFillLast(bs.y0(), bs.y0Error(), fitNLumi_);
1515     hs[k_z0_lumi]->ShiftFillLast(bs.z0(), bs.z0Error(), fitNLumi_);
1516   }
1517 
1518   // Fill summary report
1519   if (countFitting) {
1520     for (int n = 0; n < nFitElements_; n++) {
1521       reportSummaryContents[n]->Fill(summaryContent_[n] / (float)nFits_);
1522     }
1523 
1524     summarySum_ = 0;
1525     for (int ii = 0; ii < nFitElements_; ii++) {
1526       summarySum_ += summaryContent_[ii];
1527     }
1528     reportSummary_ = summarySum_ / (nFitElements_ * nFits_);
1529     if (reportSummary)
1530       reportSummary->Fill(reportSummary_);
1531 
1532     for (int bi = 0; bi < nFitElements_; bi++) {
1533       reportSummaryMap->setBinContent(1, bi + 1, summaryContent_[bi] / (float)nFits_);
1534     }
1535   }
1536 
1537   if ((resetFitNLumi_ > 0 &&
1538        ((onlineMode_ &&
1539          countLumi_ == resetFitNLumi_) ||  //OR it should be currentLumi_ (if in sequence then does not mattar)
1540         (!onlineMode_ && countLumi_ == resetFitNLumi_))) ||
1541       (StartAverage_)) {
1542     edm::LogInfo("BeamMonitor") << "FitAndFill:: The flag is ON for running average Beam Spot fit" << endl;
1543     StartAverage_ = true;
1544     firstAverageFit_++;
1545     resetHistos_ = true;
1546     nthBSTrk_ = 0;
1547     beginLumiOfBSFit_ = 0;
1548     refBStime[0] = 0;
1549   }
1550 }
1551 
1552 //--------------------------------------------------------
1553 void BeamMonitor::RestartFitting() {
1554   if (debug_)
1555     edm::LogInfo("BeamMonitor")
1556         << " RestartingFitting:: Restart Beami everything to a fresh start !!! because Gap is > 10 LS" << endl;
1557   //track based fit reset here
1558   resetHistos_ = true;
1559   nthBSTrk_ = 0;
1560   theBeamFitter->resetTrkVector();
1561   theBeamFitter->resetLSRange();
1562   theBeamFitter->resetRefTime();
1563   theBeamFitter->resetPVFitter();
1564   theBeamFitter->resetCutFlow();
1565   beginLumiOfBSFit_ = 0;
1566   refBStime[0] = 0;
1567   //pv based fit iis reset here
1568   h_PVx[0]->Reset();
1569   h_PVy[0]->Reset();
1570   h_PVz[0]->Reset();
1571   beginLumiOfPVFit_ = 0;
1572   refPVtime[0] = 0;
1573   //Clear all the Maps here
1574   mapPVx.clear();
1575   mapPVy.clear();
1576   mapPVz.clear();
1577   mapNPV.clear();
1578   mapBeginBSLS.clear();
1579   mapBeginPVLS.clear();
1580   mapBeginBSTime.clear();
1581   mapBeginPVTime.clear();
1582   mapLSBSTrkSize.clear();
1583   mapLSPVStoreSize.clear();
1584   mapLSCF.clear();
1585   countGapLumi_ = 0;
1586   countLumi_ = 0;
1587   StartAverage_ = false;
1588 }
1589 
1590 //-------------------------------------------------------
1591 void BeamMonitor::dqmEndRun(const Run& r, const EventSetup& context) {
1592   if (debug_)
1593     edm::LogInfo("BeamMonitor") << "dqmEndRun:: Clearing all the Maps " << endl;
1594   //Clear all the Maps here
1595   mapPVx.clear();
1596   mapPVy.clear();
1597   mapPVz.clear();
1598   mapNPV.clear();
1599   mapBeginBSLS.clear();
1600   mapBeginPVLS.clear();
1601   mapBeginBSTime.clear();
1602   mapBeginPVTime.clear();
1603   mapLSBSTrkSize.clear();
1604   mapLSPVStoreSize.clear();
1605   mapLSCF.clear();
1606 
1607   if (useLockRecords_ && onlineDbService_.isAvailable()) {
1608     onlineDbService_->releaseLocks();
1609   }
1610 }
1611 //--------------------------------------------------------
1612 void BeamMonitor::scrollTH1(TH1* h, time_t ref) {
1613   char offsetTime[64];
1614   formatFitTime(offsetTime, ref);
1615   TDatime da(offsetTime);
1616   if (lastNZbin > 0) {
1617     double val = h->GetBinContent(lastNZbin);
1618     double valErr = h->GetBinError(lastNZbin);
1619     h->Reset();
1620     h->GetXaxis()->SetTimeOffset(da.Convert(kTRUE));
1621     int bin = (lastNZbin > buffTime ? buffTime : 1);
1622     h->SetBinContent(bin, val);
1623     h->SetBinError(bin, valErr);
1624   } else {
1625     h->Reset();
1626     h->GetXaxis()->SetTimeOffset(da.Convert(kTRUE));
1627   }
1628 }
1629 
1630 //--------------------------------------------------------
1631 // Method to check whether to chane histogram time offset (forward only)
1632 bool BeamMonitor::testScroll(time_t& tmpTime_, time_t& refTime_) {
1633   bool scroll_ = false;
1634   if (tmpTime_ - refTime_ >= intervalInSec_) {
1635     scroll_ = true;
1636     edm::LogInfo("BeamMonitor") << "testScroll::  Reset Time Offset" << std::endl;
1637     lastNZbin = intervalInSec_;
1638     for (int bin = intervalInSec_; bin >= 1; bin--) {
1639       if (hs[k_x0_time]->getBinContent(bin) > 0) {
1640         lastNZbin = bin;
1641         break;
1642       }
1643     }
1644     edm::LogInfo("BeamMonitor") << "testScroll::  Last non zero bin = " << lastNZbin << std::endl;
1645     if (tmpTime_ - refTime_ >= intervalInSec_ + lastNZbin) {
1646       edm::LogInfo("BeamMonitor") << "testScroll::  Time difference too large since last readout" << std::endl;
1647       lastNZbin = 0;
1648       refTime_ = tmpTime_ - buffTime;
1649     } else {
1650       edm::LogInfo("BeamMonitor") << "testScroll::  Offset to last record" << std::endl;
1651       int offset = ((lastNZbin > buffTime) ? (lastNZbin - buffTime) : (lastNZbin - 1));
1652       refTime_ += offset;
1653     }
1654   }
1655   return scroll_;
1656 }
1657 
1658 DEFINE_FWK_MODULE(BeamMonitor);
1659 
1660 // Local Variables:
1661 // show-trailing-whitespace: t
1662 // truncate-lines: t
1663 // End:
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!