From 98288b8123bddc53db62035ddd995078eece82a8 Mon Sep 17 00:00:00 2001 From: Mauro Date: Thu, 24 Sep 2015 00:38:47 +0200 Subject: [PATCH] Now the code is able to handle non consecutive LS --- DQM/BeamMonitor/plugins/Vx3DHLTAnalyzer.cc | 821 ++++++++++-------- DQM/BeamMonitor/plugins/Vx3DHLTAnalyzer.h | 231 +++-- .../beampixel_dqm_sourceclient-live_cfg.py | 112 +-- 3 files changed, 619 insertions(+), 545 deletions(-) diff --git a/DQM/BeamMonitor/plugins/Vx3DHLTAnalyzer.cc b/DQM/BeamMonitor/plugins/Vx3DHLTAnalyzer.cc index 441c391e0e649..4a9d721acbb88 100644 --- a/DQM/BeamMonitor/plugins/Vx3DHLTAnalyzer.cc +++ b/DQM/BeamMonitor/plugins/Vx3DHLTAnalyzer.cc @@ -1,16 +1,16 @@ // -*- C++ -*- -// -// Package: Vx3DHLTAnalyzer -// Class: Vx3DHLTAnalyzer -// -/**\class Vx3DHLTAnalyzer Vx3DHLTAnalyzer.cc plugins/Vx3DHLTAnalyzer.cc - - Description: beam-spot monitor entirely based on pixel detector information - Implementation: the monitoring is based on a 3D fit to the vertex cloud +// Package: Vx3DHLTAnalyzer +// Class: Vx3DHLTAnalyzer + +/* + Class Vx3DHLTAnalyzer Vx3DHLTAnalyzer.cc plugins/Vx3DHLTAnalyzer.cc + + Description: beam-spot monitor entirely based on pixel detector information + Implementation: the monitoring is based on a 3D fit to the vertex cloud */ -// -// Original Author: Mauro Dinardo, 28 S-012, +41-22-767-8302, -// Created: Tue Feb 23 13:15:31 CET 2010 + +// Original Author: Mauro Dinardo, 28 S-012, +41-22-767-8302 +// Created: Tue Feb 23 13:15:31 CET 2010 #include "DQM/BeamMonitor/plugins/Vx3DHLTAnalyzer.h" @@ -23,26 +23,40 @@ #include +// ### Calling namespaces ### +using namespace std; +using namespace edm; +using namespace reco; + + Vx3DHLTAnalyzer::Vx3DHLTAnalyzer (const ParameterSet& iConfig) { - debugMode = true; - nLumiReset = 2; // Number of integrated lumis to perform the fit - dataFromFit = true; // The Beam Spot data can be either taken from the histograms or from the fit results - minNentries = 20; // Minimum number of good vertices to perform the fit - xRange = 1.; // [cm] - xStep = 0.001; // [cm] - yRange = 1.; // [cm] - yStep = 0.001; // [cm] - zRange = 30.; // [cm] - zStep = 0.05; // [cm] - VxErrCorr = 1.3; - fileName = "BeamPixelResults.txt"; + debugMode = true; + nLumiFit = 2; // Number of integrated lumis to perform the fit + maxLumiIntegration = 15; // If failing fits, this is the maximum number of integrated lumis after which a reset is issued + nLumiXaxisRange = 3000; // Correspond to about 20h of data taking: 20h * 60min * 60s / 23s per lumi-block = 3130 + dataFromFit = true; // The Beam Spot data can be either taken from the histograms or from the fit results + minNentries = 20; // Minimum number of good vertices to perform the fit + xRange = 0.8; // [cm] + xStep = 0.001; // [cm] + yRange = 0.8; // [cm] + yStep = 0.001; // [cm] + zRange = 30.; // [cm] + zStep = 0.04; // [cm] + VxErrCorr = 1.3; + minVxDoF = 10.; // Good-vertex selection cut + // For vertex fitter without track-weight: d.o.f. = 2*NTracks - 3 + // For vertex fitter with track-weight: d.o.f. = sum_NTracks(2*track_weight) - 3 + minVxWgt = 0.5; // Good-vertex selection cut + fileName = "BeamPixelResults.txt"; vertexCollection = consumes (iConfig.getUntrackedParameter("vertexCollection", InputTag("pixelVertices"))); pixelHitCollection = consumes(iConfig.getUntrackedParameter("pixelHitCollection", InputTag("siPixelRecHits"))); debugMode = iConfig.getParameter("debugMode"); - nLumiReset = iConfig.getParameter("nLumiReset"); + nLumiFit = iConfig.getParameter("nLumiFit"); + maxLumiIntegration = iConfig.getParameter("maxLumiIntegration"); + nLumiXaxisRange = iConfig.getParameter("nLumiXaxisRange"); dataFromFit = iConfig.getParameter("dataFromFit"); minNentries = iConfig.getParameter("minNentries"); xRange = iConfig.getParameter("xRange"); @@ -52,12 +66,23 @@ Vx3DHLTAnalyzer::Vx3DHLTAnalyzer (const ParameterSet& iConfig) zRange = iConfig.getParameter("zRange"); zStep = iConfig.getParameter("zStep"); VxErrCorr = iConfig.getParameter("VxErrCorr"); + minVxDoF = iConfig.getParameter("minVxDoF"); + minVxWgt = iConfig.getParameter("minVxWgt"); fileName = iConfig.getParameter("fileName"); + + + // ### Set internal variables ### + nParams = 9; // Number of free parameters in the fit + internalDebug = false; + considerVxCovariance = true; // Deconvolute vertex covariance matrix + pi = 3.141592653589793238; + // ############################## } Vx3DHLTAnalyzer::~Vx3DHLTAnalyzer () { + reset("scratch"); } @@ -87,18 +112,35 @@ void Vx3DHLTAnalyzer::analyze (const Event& iEvent, const EventSetup& iSetup) outputDebugFile.close(); outputDebugFile.open(debugFile.str().c_str(), ios::app); } - + beginLuminosityBlock(iEvent.getLuminosityBlock(),iSetup); } else if (beginTimeOfFit != 0) { totalHits += HitCounter(iEvent); + if (internalDebug == true) + { + cout << "[Vx3DHLTAnalyzer]::\tI found " << totalHits << " pixel hits until now" << endl; + cout << "[Vx3DHLTAnalyzer]::\tIn this event there are " << Vx3DCollection->size() << " vertex cadidates" << endl; + } + for (vector::const_iterator it3DVx = Vx3DCollection->begin(); it3DVx != Vx3DCollection->end(); it3DVx++) { - if ((it3DVx->isValid() == true) && - (it3DVx->isFake() == false) && - (it3DVx->ndof() >= minVxDoF)) + if (internalDebug == true) + { + cout << "[Vx3DHLTAnalyzer]::\tVertex selections:" << endl; + cout << "[Vx3DHLTAnalyzer]::\tisValid = " << it3DVx->isValid() << endl; + cout << "[Vx3DHLTAnalyzer]::\tisFake = " << it3DVx->isFake() << endl; + cout << "[Vx3DHLTAnalyzer]::\tnodof = " << it3DVx->ndof() << endl; + cout << "[Vx3DHLTAnalyzer]::\ttracksSize = " << it3DVx->tracksSize() << endl; + } + + if ((it3DVx->isValid() == true) && + (it3DVx->isFake() == false) && + (it3DVx->ndof() >= minVxDoF) && + (it3DVx->tracksSize() > 0) && + ((it3DVx->ndof()+3.) / ((double)it3DVx->tracksSize()) >= 2.*minVxWgt)) { for (i = 0; i < DIM; i++) { @@ -107,34 +149,42 @@ void Vx3DHLTAnalyzer::analyze (const Event& iEvent, const EventSetup& iSetup) MyVertex.Covariance[i][j] = it3DVx->covariance(i,j); if (isNotFinite(MyVertex.Covariance[i][j]) == true) break; } + if (j != DIM) break; } - det = std::fabs(MyVertex.Covariance[0][0])*(std::fabs(MyVertex.Covariance[1][1])*std::fabs(MyVertex.Covariance[2][2]) - MyVertex.Covariance[1][2]*MyVertex.Covariance[1][2]) - - MyVertex.Covariance[0][1]*(MyVertex.Covariance[0][1]*std::fabs(MyVertex.Covariance[2][2]) - MyVertex.Covariance[0][2]*MyVertex.Covariance[1][2]) + - MyVertex.Covariance[0][2]*(MyVertex.Covariance[0][1]*MyVertex.Covariance[1][2] - MyVertex.Covariance[0][2]*std::fabs(MyVertex.Covariance[1][1])); + + if (i == DIM) + det = std::fabs(MyVertex.Covariance[0][0])*(std::fabs(MyVertex.Covariance[1][1])*std::fabs(MyVertex.Covariance[2][2]) - MyVertex.Covariance[1][2]*MyVertex.Covariance[1][2]) - + MyVertex.Covariance[0][1]*(MyVertex.Covariance[0][1]*std::fabs(MyVertex.Covariance[2][2]) - MyVertex.Covariance[0][2]*MyVertex.Covariance[1][2]) + + MyVertex.Covariance[0][2]*(MyVertex.Covariance[0][1]*MyVertex.Covariance[1][2] - MyVertex.Covariance[0][2]*std::fabs(MyVertex.Covariance[1][1])); + if ((i == DIM) && (det > 0.)) { + if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tVertex accepted !" << endl; + MyVertex.x = it3DVx->x(); MyVertex.y = it3DVx->y(); MyVertex.z = it3DVx->z(); Vertices.push_back(MyVertex); + + Vx_X->Fill(it3DVx->x()); + Vx_Y->Fill(it3DVx->y()); + Vx_Z->Fill(it3DVx->z()); + + Vx_ZX->Fill(it3DVx->z(), it3DVx->x()); + Vx_ZY->Fill(it3DVx->z(), it3DVx->y()); + Vx_XY->Fill(it3DVx->x(), it3DVx->y()); } else if (internalDebug == true) { - cout << "Vertex discarded !" << endl; + cout << "[Vx3DHLTAnalyzer]::\tVertex discarded !" << endl; + for (i = 0; i < DIM; i++) for (j = 0; j < DIM; j++) cout << "(i,j) --> " << i << "," << j << " --> " << MyVertex.Covariance[i][j] << endl; } - - Vx_X->Fill(it3DVx->x()); - Vx_Y->Fill(it3DVx->y()); - Vx_Z->Fill(it3DVx->z()); - - Vx_ZX->Fill(it3DVx->z(), it3DVx->x()); - Vx_ZY->Fill(it3DVx->z(), it3DVx->y()); - Vx_XY->Fill(it3DVx->x(), it3DVx->y()); } + else if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tVertex discarded !" << endl; } } } @@ -165,7 +215,7 @@ string Vx3DHLTAnalyzer::formatTime (const time_t& t) } -double Gauss3DFunc(const double* par) +double Gauss3DFunc (const double* par) { double K[DIM][DIM]; // Covariance Matrix double M[DIM][DIM]; // K^-1 @@ -236,22 +286,26 @@ double Gauss3DFunc(const double* par) int Vx3DHLTAnalyzer::MyFit (vector* vals) { - // RETURN CODE: - // 0 == OK - // -2 == NO OK - not enough "minNentries" - // Any other number == NO OK - unsigned int nParams = 9; - + // ############################################ + // # RETURN CODE: # + // # >0 == NO OK - fit status (MINUIT manual) # + // # 0 == OK # + // # -1 == NO OK - not finite edm # + // # -2 == NO OK - not enough "minNentries" # + // # -3 == NO OK - not finite errors # + // # -4 == NO OK - negative determinant # + // # -5 == NO OK - maxLumiIntegration reached # + // ############################################ + if ((vals != NULL) && (vals->size() == nParams*2)) { - double nSigmaXY = 100.; - double nSigmaZ = 100.; - double varFactor = 4./25.; // Take into account the difference between the RMS and sigma (RMS usually greater than sigma) - double parDistanceXY = 0.005; // Unit: [cm] - double parDistanceZ = 0.5; // Unit: [cm] - double parDistanceddZ = 1e-3; // Unit: [rad] - double parDistanceCxy = 1e-5; // Unit: [cm^2] - double bestEdm = 1e-1; + double nSigmaXY = 10.; + double nSigmaZ = 10.; + double parDistanceXY = 1e-3; // Unit: [cm] + double parDistanceZ = 1e-2; // Unit: [cm] + double parDistanceddZ = 1e-3; // Unit: [rad] + double parDistanceCxy = 1e-5; // Unit: [cm^2] + double bestEdm; const unsigned int trials = 4; double largerDist[trials] = {0.1, 5., 10., 100.}; @@ -268,28 +322,27 @@ int Vx3DHLTAnalyzer::MyFit (vector* vals) vector::const_iterator it = vals->begin(); ROOT::Math::Minimizer* Gauss3D = ROOT::Math::Factory::CreateMinimizer("Minuit2","Migrad"); - Gauss3D->SetMaxFunctionCalls(1e4); - Gauss3D->SetTolerance(1e-9); // Tolerance on likelihood + Gauss3D->SetErrorDef(1.0); if (internalDebug == true) Gauss3D->SetPrintLevel(3); else Gauss3D->SetPrintLevel(0); ROOT::Math::Functor _Gauss3DFunc(&Gauss3DFunc,nParams); Gauss3D->SetFunction(_Gauss3DFunc); - if (internalDebug == true) cout << "\n@@@ START FITTING @@@" << endl; + if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\t@@@ START FITTING @@@" << endl; // @@@ Fit at X-deltaMean | X | X+deltaMean @@@ bestEdm = 1.; for (int i = 0; i < 3; i++) { - deltaMean = (double(i)-1.)*std::sqrt((*(it+0))*varFactor); - if (internalDebug == true) cout << "deltaMean --> " << deltaMean << endl; + deltaMean = (double(i)-1.)*std::sqrt(*(it+0)); + if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tdeltaMean --> " << deltaMean << endl; Gauss3D->Clear(); - Gauss3D->SetVariable(0,"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY); - Gauss3D->SetVariable(1,"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY); - Gauss3D->SetVariable(2,"var z ", *(it+2), parDistanceZ*parDistanceZ); + Gauss3D->SetVariable(0,"var x ", *(it+0), parDistanceXY * parDistanceXY); + Gauss3D->SetVariable(1,"var y ", *(it+1), parDistanceXY * parDistanceXY); + Gauss3D->SetVariable(2,"var z ", *(it+2), parDistanceZ * parDistanceZ); Gauss3D->SetVariable(3,"cov xy", *(it+3), parDistanceCxy); Gauss3D->SetVariable(4,"dydz ", *(it+4), parDistanceddZ); Gauss3D->SetVariable(5,"dxdz ", *(it+5), parDistanceddZ); @@ -311,8 +364,14 @@ int Vx3DHLTAnalyzer::MyFit (vector* vals) edm = Gauss3D->Edm(); if (counterVx < minNentries) goodData = -2; - else if (isNotFinite(edm) == true) goodData = -1; - else for (unsigned int j = 0; j < nParams; j++) if (isNotFinite(Gauss3D->Errors()[j]) == true) { goodData = -1; break; } + else if (isNotFinite(edm) == true) { goodData = -1; if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl; } + else for (unsigned int j = 0; j < nParams; j++) + if (isNotFinite(Gauss3D->Errors()[j]) == true) + { + goodData = -3; + if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl; + break; + } if (goodData == 0) { covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3]; @@ -321,33 +380,33 @@ int Vx3DHLTAnalyzer::MyFit (vector* vals) det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) - Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) + covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1])); - if (det < 0.) { goodData = -1; if (internalDebug == true) cout << "Negative determinant !" << endl; } + if (det < 0.) { goodData = -4; if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl; } } if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementX = i; } } - if (internalDebug == true) cout << "Found bestMovementX --> " << bestMovementX << endl; + if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tFound bestMovementX --> " << bestMovementX << endl; // @@@ Fit at Y-deltaMean | Y | Y+deltaMean @@@ bestEdm = 1.; for (int i = 0; i < 3; i++) { - deltaMean = (double(i)-1.)*std::sqrt((*(it+1))*varFactor); + deltaMean = (double(i)-1.)*std::sqrt(*(it+1)); if (internalDebug == true) { - cout << "deltaMean --> " << deltaMean << endl; - cout << "deltaMean X --> " << (double(bestMovementX)-1.)*std::sqrt((*(it+0))*varFactor) << endl; + cout << "[Vx3DHLTAnalyzer]::\tdeltaMean --> " << deltaMean << endl; + cout << "[Vx3DHLTAnalyzer]::\tdeltaMean X --> " << (double(bestMovementX)-1.)*std::sqrt(*(it+0)) << endl; } Gauss3D->Clear(); - Gauss3D->SetVariable(0,"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY); - Gauss3D->SetVariable(1,"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY); - Gauss3D->SetVariable(2,"var z ", *(it+2), parDistanceZ*parDistanceZ); + Gauss3D->SetVariable(0,"var x ", *(it+0), parDistanceXY * parDistanceXY); + Gauss3D->SetVariable(1,"var y ", *(it+1), parDistanceXY * parDistanceXY); + Gauss3D->SetVariable(2,"var z ", *(it+2), parDistanceZ * parDistanceZ); Gauss3D->SetVariable(3,"cov xy", *(it+3), parDistanceCxy); Gauss3D->SetVariable(4,"dydz ", *(it+4), parDistanceddZ); Gauss3D->SetVariable(5,"dxdz ", *(it+5), parDistanceddZ); - Gauss3D->SetVariable(6,"mean x", *(it+6)+(double(bestMovementX)-1.)*std::sqrt((*(it+0))*varFactor), parDistanceXY); + Gauss3D->SetVariable(6,"mean x", *(it+6)+(double(bestMovementX)-1.)*std::sqrt(*(it+0)), parDistanceXY); Gauss3D->SetVariable(7,"mean y", *(it+7)+deltaMean, parDistanceXY); Gauss3D->SetVariable(8,"mean z", *(it+8), parDistanceZ); @@ -365,8 +424,14 @@ int Vx3DHLTAnalyzer::MyFit (vector* vals) edm = Gauss3D->Edm(); if (counterVx < minNentries) goodData = -2; - else if (isNotFinite(edm) == true) goodData = -1; - else for (unsigned int j = 0; j < nParams; j++) if (isNotFinite(Gauss3D->Errors()[j]) == true) { goodData = -1; break; } + else if (isNotFinite(edm) == true) { goodData = -1; if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl; } + else for (unsigned int j = 0; j < nParams; j++) + if (isNotFinite(Gauss3D->Errors()[j]) == true) + { + goodData = -3; + if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl; + break; + } if (goodData == 0) { covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3]; @@ -375,12 +440,12 @@ int Vx3DHLTAnalyzer::MyFit (vector* vals) det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) - Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) + covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1])); - if (det < 0.) { goodData = -1; if (internalDebug == true) cout << "Negative determinant !" << endl; } + if (det < 0.) { goodData = -4; if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl; } } if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementY = i; } } - if (internalDebug == true) cout << "Found bestMovementY --> " << bestMovementY << endl; + if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tFound bestMovementY --> " << bestMovementY << endl; // @@@ Fit at Z-deltaMean | Z | Z+deltaMean @@@ bestEdm = 1.; @@ -389,21 +454,21 @@ int Vx3DHLTAnalyzer::MyFit (vector* vals) deltaMean = (double(i)-1.)*std::sqrt(*(it+2)); if (internalDebug == true) { - cout << "deltaMean --> " << deltaMean << endl; - cout << "deltaMean X --> " << (double(bestMovementX)-1.)*std::sqrt((*(it+0))*varFactor) << endl; - cout << "deltaMean Y --> " << (double(bestMovementY)-1.)*std::sqrt((*(it+1))*varFactor) << endl; + cout << "[Vx3DHLTAnalyzer]::\tdeltaMean --> " << deltaMean << endl; + cout << "[Vx3DHLTAnalyzer]::\tdeltaMean X --> " << (double(bestMovementX)-1.)*std::sqrt(*(it+0)) << endl; + cout << "[Vx3DHLTAnalyzer]::\tdeltaMean Y --> " << (double(bestMovementY)-1.)*std::sqrt(*(it+1)) << endl; } Gauss3D->Clear(); - Gauss3D->SetVariable(0,"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY); - Gauss3D->SetVariable(1,"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY); - Gauss3D->SetVariable(2,"var z ", *(it+2), parDistanceZ*parDistanceZ); + Gauss3D->SetVariable(0,"var x ", *(it+0), parDistanceXY * parDistanceXY); + Gauss3D->SetVariable(1,"var y ", *(it+1), parDistanceXY * parDistanceXY); + Gauss3D->SetVariable(2,"var z ", *(it+2), parDistanceZ * parDistanceZ); Gauss3D->SetVariable(3,"cov xy", *(it+3), parDistanceCxy); Gauss3D->SetVariable(4,"dydz ", *(it+4), parDistanceddZ); Gauss3D->SetVariable(5,"dxdz ", *(it+5), parDistanceddZ); - Gauss3D->SetVariable(6,"mean x", *(it+6)+(double(bestMovementX)-1.)*std::sqrt((*(it+0))*varFactor), parDistanceXY); - Gauss3D->SetVariable(7,"mean y", *(it+7)+(double(bestMovementY)-1.)*std::sqrt((*(it+1))*varFactor), parDistanceXY); + Gauss3D->SetVariable(6,"mean x", *(it+6)+(double(bestMovementX)-1.)*std::sqrt(*(it+0)), parDistanceXY); + Gauss3D->SetVariable(7,"mean y", *(it+7)+(double(bestMovementY)-1.)*std::sqrt(*(it+1)), parDistanceXY); Gauss3D->SetVariable(8,"mean z", *(it+8)+deltaMean, parDistanceZ); // Set the central positions of the centroid for vertex rejection @@ -420,8 +485,14 @@ int Vx3DHLTAnalyzer::MyFit (vector* vals) edm = Gauss3D->Edm(); if (counterVx < minNentries) goodData = -2; - else if (isNotFinite(edm) == true) goodData = -1; - else for (unsigned int j = 0; j < nParams; j++) if (isNotFinite(Gauss3D->Errors()[j]) == true) { goodData = -1; break; } + else if (isNotFinite(edm) == true) { goodData = -1; if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl; } + else for (unsigned int j = 0; j < nParams; j++) + if (isNotFinite(Gauss3D->Errors()[j]) == true) + { + goodData = -3; + if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl; + break; + } if (goodData == 0) { covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3]; @@ -430,24 +501,24 @@ int Vx3DHLTAnalyzer::MyFit (vector* vals) det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) - Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) + covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1])); - if (det < 0.) { goodData = -1; if (internalDebug == true) cout << "Negative determinant !" << endl; } + if (det < 0.) { goodData = -4; if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl; } } if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementZ = i; } } - if (internalDebug == true) cout << "Found bestMovementZ --> " << bestMovementZ << endl; + if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tFound bestMovementZ --> " << bestMovementZ << endl; Gauss3D->Clear(); // @@@ FINAL FIT @@@ - Gauss3D->SetVariable(0,"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY); - Gauss3D->SetVariable(1,"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY); - Gauss3D->SetVariable(2,"var z ", *(it+2), parDistanceZ*parDistanceZ); + Gauss3D->SetVariable(0,"var x ", *(it+0), parDistanceXY * parDistanceXY); + Gauss3D->SetVariable(1,"var y ", *(it+1), parDistanceXY * parDistanceXY); + Gauss3D->SetVariable(2,"var z ", *(it+2), parDistanceZ * parDistanceZ); Gauss3D->SetVariable(3,"cov xy", *(it+3), parDistanceCxy); Gauss3D->SetVariable(4,"dydz ", *(it+4), parDistanceddZ); Gauss3D->SetVariable(5,"dxdz ", *(it+5), parDistanceddZ); - Gauss3D->SetVariable(6,"mean x", *(it+6)+(double(bestMovementX)-1.)*std::sqrt((*(it+0))*varFactor), parDistanceXY); - Gauss3D->SetVariable(7,"mean y", *(it+7)+(double(bestMovementY)-1.)*std::sqrt((*(it+1))*varFactor), parDistanceXY); + Gauss3D->SetVariable(6,"mean x", *(it+6)+(double(bestMovementX)-1.)*std::sqrt(*(it+0)), parDistanceXY); + Gauss3D->SetVariable(7,"mean y", *(it+7)+(double(bestMovementY)-1.)*std::sqrt(*(it+1)), parDistanceXY); Gauss3D->SetVariable(8,"mean z", *(it+8)+(double(bestMovementZ)-1.)*std::sqrt(*(it+2)), parDistanceZ); // Set the central positions of the centroid for vertex rejection @@ -464,8 +535,14 @@ int Vx3DHLTAnalyzer::MyFit (vector* vals) edm = Gauss3D->Edm(); if (counterVx < minNentries) goodData = -2; - else if (isNotFinite(edm) == true) goodData = -1; - else for (unsigned int j = 0; j < nParams; j++) if (isNotFinite(Gauss3D->Errors()[j]) == true) { goodData = -1; break; } + else if (isNotFinite(edm) == true) { goodData = -1; if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl; } + else for (unsigned int j = 0; j < nParams; j++) + if (isNotFinite(Gauss3D->Errors()[j]) == true) + { + goodData = -3; + if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl; + break; + } if (goodData == 0) { covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3]; @@ -474,7 +551,7 @@ int Vx3DHLTAnalyzer::MyFit (vector* vals) det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) - Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) + covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1])); - if (det < 0.) { goodData = -1; if (internalDebug == true) cout << "Negative determinant !" << endl; } + if (det < 0.) { goodData = -4; if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl; } } // @@@ FIT WITH DIFFERENT PARAMETER DISTANCES @@@ @@ -484,17 +561,17 @@ int Vx3DHLTAnalyzer::MyFit (vector* vals) { Gauss3D->Clear(); - if (internalDebug == true) cout << "FIT WITH DIFFERENT PARAMETER DISTANCES - STEP " << i+1 << endl; + if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tFIT WITH DIFFERENT PARAMETER DISTANCES - STEP " << i+1 << endl; - Gauss3D->SetVariable(0,"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY * largerDist[i]); - Gauss3D->SetVariable(1,"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY * largerDist[i]); - Gauss3D->SetVariable(2,"var z ", *(it+2), parDistanceZ*parDistanceZ * largerDist[i]); + Gauss3D->SetVariable(0,"var x ", *(it+0), parDistanceXY * parDistanceXY * largerDist[i]); + Gauss3D->SetVariable(1,"var y ", *(it+1), parDistanceXY * parDistanceXY * largerDist[i]); + Gauss3D->SetVariable(2,"var z ", *(it+2), parDistanceZ * parDistanceZ * largerDist[i]); Gauss3D->SetVariable(3,"cov xy", *(it+3), parDistanceCxy * largerDist[i]); Gauss3D->SetVariable(4,"dydz ", *(it+4), parDistanceddZ * largerDist[i]); Gauss3D->SetVariable(5,"dxdz ", *(it+5), parDistanceddZ * largerDist[i]); - Gauss3D->SetVariable(6,"mean x", *(it+6)+(double(bestMovementX)-1.)*std::sqrt((*(it+0))*varFactor), parDistanceXY * largerDist[i]); - Gauss3D->SetVariable(7,"mean y", *(it+7)+(double(bestMovementY)-1.)*std::sqrt((*(it+1))*varFactor), parDistanceXY * largerDist[i]); - Gauss3D->SetVariable(8,"mean z", *(it+8)+(double(bestMovementZ)-1.)*std::sqrt(*(it+2)), parDistanceZ * largerDist[i]); + Gauss3D->SetVariable(6,"mean x", *(it+6)+(double(bestMovementX)-1.)*std::sqrt(*(it+0)), parDistanceXY * largerDist[i]); + Gauss3D->SetVariable(7,"mean y", *(it+7)+(double(bestMovementY)-1.)*std::sqrt(*(it+1)), parDistanceXY * largerDist[i]); + Gauss3D->SetVariable(8,"mean z", *(it+8)+(double(bestMovementZ)-1.)*std::sqrt(*(it+2)), parDistanceZ * largerDist[i]); // Set the central positions of the centroid for vertex rejection xPos = Gauss3D->X()[6]; @@ -510,8 +587,14 @@ int Vx3DHLTAnalyzer::MyFit (vector* vals) edm = Gauss3D->Edm(); if (counterVx < minNentries) goodData = -2; - else if (isNotFinite(edm) == true) goodData = -1; - else for (unsigned int j = 0; j < nParams; j++) if (isNotFinite(Gauss3D->Errors()[j]) == true) { goodData = -1; break; } + else if (isNotFinite(edm) == true) { goodData = -1; if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl; } + else for (unsigned int j = 0; j < nParams; j++) + if (isNotFinite(Gauss3D->Errors()[j]) == true) + { + goodData = -3; + if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl; + break; + } if (goodData == 0) { covyz = Gauss3D->X()[4]*(std::fabs(Gauss3D->X()[2])-std::fabs(Gauss3D->X()[1])) - Gauss3D->X()[5]*Gauss3D->X()[3]; @@ -520,7 +603,7 @@ int Vx3DHLTAnalyzer::MyFit (vector* vals) det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1])*std::fabs(Gauss3D->X()[2]) - covyz*covyz) - Gauss3D->X()[3] * (Gauss3D->X()[3]*std::fabs(Gauss3D->X()[2]) - covxz*covyz) + covxz * (Gauss3D->X()[3]*covyz - covxz*std::fabs(Gauss3D->X()[1])); - if (det < 0.) { goodData = -1; if (internalDebug == true) cout << "Negative determinant !" << endl; } + if (det < 0.) { goodData = -4; if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl; } } } else break; } @@ -542,6 +625,18 @@ int Vx3DHLTAnalyzer::MyFit (vector* vals) void Vx3DHLTAnalyzer::reset (string ResetType) { + if ((debugMode == true) && (outputDebugFile.is_open() == true)) + { + outputDebugFile << "Runnumber " << runNumber << endl; + outputDebugFile << "BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) << " " << (beginTimeOfFit >> 32) << endl; + outputDebugFile << "BeginLumiRange " << beginLumiOfFit << endl; + outputDebugFile << "EndTimeOfFit " << formatTime(endTimeOfFit >> 32) << " " << (endTimeOfFit >> 32) << endl; + outputDebugFile << "EndLumiRange " << endLumiOfFit << endl; + outputDebugFile << "LumiCounter " << lumiCounter << endl; + outputDebugFile << "LastLumiOfFit " << lastLumiOfFit << endl; + } + + if (ResetType.compare("scratch") == 0) { runNumber = 0; @@ -569,73 +664,54 @@ void Vx3DHLTAnalyzer::reset (string ResetType) dydzlumi->Reset(); hitCounter->Reset(); - hitCountHistory->Reset(); goodVxCounter->Reset(); - goodVxCountHistory->Reset(); + statusCounter->Reset(); fitResults->Reset(); - reportSummary->Fill(-1.); - reportSummaryMap->Fill(0.5, 0.5, -1.); + reportSummary->Fill(-1); + reportSummaryMap->getTH1()->SetBinContent(1, 1, -1); Vertices.clear(); - lumiCounter = 0; - lumiCounterHisto = 0; - totalHits = 0; - beginTimeOfFit = 0; - endTimeOfFit = 0; - beginLumiOfFit = 0; - endLumiOfFit = 0; + lumiCounter = 0; + totalHits = 0; + beginTimeOfFit = 0; + endTimeOfFit = 0; + beginLumiOfFit = 0; + endLumiOfFit = 0; + + if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tReset issued: scratch" << endl; + if ((debugMode == true) && (outputDebugFile.is_open() == true)) outputDebugFile << "Reset -scratch- issued\n" << endl; } else if (ResetType.compare("whole") == 0) { Vx_X->Reset(); Vx_Y->Reset(); Vx_Z->Reset(); - + Vx_ZX->Reset(); Vx_ZY->Reset(); Vx_XY->Reset(); Vertices.clear(); - lumiCounter = 0; - lumiCounterHisto = 0; - totalHits = 0; - beginTimeOfFit = 0; - endTimeOfFit = 0; - beginLumiOfFit = 0; - endLumiOfFit = 0; + lumiCounter = 0; + totalHits = 0; + beginTimeOfFit = 0; + endTimeOfFit = 0; + beginLumiOfFit = 0; + endLumiOfFit = 0; + + if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tReset issued: whole" << endl; + if ((debugMode == true) && (outputDebugFile.is_open() == true)) outputDebugFile << "Reset -whole- issued\n" << endl; } - else if (ResetType.compare("partial") == 0) - { - Vx_X->Reset(); - Vx_Y->Reset(); - Vx_Z->Reset(); - - Vertices.clear(); - - lumiCounter = 0; - totalHits = 0; - beginTimeOfFit = 0; - endTimeOfFit = 0; - beginLumiOfFit = 0; - endLumiOfFit = 0; - } - else if (ResetType.compare("nohisto") == 0) + else if (ResetType.compare("hitCounter") == 0) { - Vertices.clear(); - - lumiCounter = 0; - lumiCounterHisto = 0; - totalHits = 0; - beginTimeOfFit = 0; - endTimeOfFit = 0; - beginLumiOfFit = 0; - endLumiOfFit = 0; + totalHits = 0; + + if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tReset issued: hitCounter" << endl; + if ((debugMode == true) && (outputDebugFile.is_open() == true)) outputDebugFile << "Reset -hitCounter- issued\n" << endl; } - else if (ResetType.compare("hitCounter") == 0) - totalHits = 0; } @@ -651,7 +727,7 @@ void Vx3DHLTAnalyzer::writeToFile (vector* vals, outputFile.open(fileName.c_str(), ios::out); - if ((outputFile.is_open() == true) && (vals != NULL) && (vals->size() == 8*2)) + if ((outputFile.is_open() == true) && (vals != NULL) && (vals->size() == (nParams-1)*2)) { vector::const_iterator it = vals->begin(); @@ -710,7 +786,7 @@ void Vx3DHLTAnalyzer::writeToFile (vector* vals, } outputFile.close(); - if ((debugMode == true) && (outputDebugFile.is_open() == true) && (vals != NULL) && (vals->size() == 8*2)) + if ((debugMode == true) && (outputDebugFile.is_open() == true) && (vals != NULL) && (vals->size() == (nParams-1)*2)) { vector::const_iterator it = vals->begin(); @@ -766,20 +842,37 @@ void Vx3DHLTAnalyzer::writeToFile (vector* vals, outputDebugFile << "EmittanceX 0.0" << endl; outputDebugFile << "EmittanceY 0.0" << endl; outputDebugFile << "BetaStar 0.0" << endl; + + outputDebugFile << "Used vertices: " << counterVx << "\n" << endl; } } +void Vx3DHLTAnalyzer::printFitParams (const vector& fitResults) +{ + cout << "var x --> " << fitResults[0] << " +/- " << fitResults[0+nParams] << endl; + cout << "var y --> " << fitResults[1] << " +/- " << fitResults[1+nParams] << endl; + cout << "var z --> " << fitResults[2] << " +/- " << fitResults[2+nParams] << endl; + cout << "cov xy --> " << fitResults[3] << " +/- " << fitResults[3+nParams] << endl; + cout << "dydz --> " << fitResults[4] << " +/- " << fitResults[4+nParams] << endl; + cout << "dxdz --> " << fitResults[5] << " +/- " << fitResults[5+nParams] << endl; + cout << "mean x --> " << fitResults[6] << " +/- " << fitResults[6+nParams] << endl; + cout << "mean y --> " << fitResults[7] << " +/- " << fitResults[7+nParams] << endl; + cout << "mean z --> " << fitResults[8] << " +/- " << fitResults[8+nParams] << endl; +} + + void Vx3DHLTAnalyzer::beginLuminosityBlock (const LuminosityBlock& lumiBlock, const EventSetup& iSetup) { + // @@@ If statement to avoid problems with non-sequential lumisections @@@ if ((lumiCounter == 0) && (lumiBlock.luminosityBlock() > lastLumiOfFit)) { beginTimeOfFit = lumiBlock.beginTime().value(); beginLumiOfFit = lumiBlock.luminosityBlock(); lumiCounter++; - lumiCounterHisto++; } - else if ((lumiCounter != 0) && (lumiBlock.luminosityBlock() >= (beginLumiOfFit+lumiCounter))) { lumiCounter++; lumiCounterHisto++; } + else if ((lumiCounter != 0) && (lumiBlock.luminosityBlock() >= (beginLumiOfFit+lumiCounter))) lumiCounter++; + else reset("scratch"); } @@ -787,22 +880,16 @@ void Vx3DHLTAnalyzer::endLuminosityBlock (const LuminosityBlock& lumiBlock, cons { stringstream histTitle; int goodData; - unsigned int nParams = 9; - if ((nLumiReset != 0) && (lumiCounter%nLumiReset == 0) && (beginTimeOfFit != 0) && (runNumber != 0)) + if ((nLumiFit != 0) && (lumiCounter%nLumiFit == 0) && (beginTimeOfFit != 0) && (runNumber != 0)) { endTimeOfFit = lumiBlock.endTime().value(); endLumiOfFit = lumiBlock.luminosityBlock(); lastLumiOfFit = endLumiOfFit; vector vals; - hitCounter->ShiftFillLast((double)totalHits, std::sqrt((double)totalHits), nLumiReset); - - if (lastLumiOfFit % prescaleHistory == 0) - { - hitCountHistory->getTH1()->SetBinContent(lastLumiOfFit, (double)totalHits); - hitCountHistory->getTH1()->SetBinError(lastLumiOfFit, std::sqrt((double)totalHits)); - } + hitCounter->getTH1()->SetBinContent(lastLumiOfFit, (double)totalHits); + hitCounter->getTH1()->SetBinError(lastLumiOfFit, std::sqrt((double)totalHits)); if (dataFromFit == true) { @@ -819,21 +906,28 @@ void Vx3DHLTAnalyzer::endLuminosityBlock (const LuminosityBlock& lumiBlock, cons fitResults.push_back(Vx_Z->getTH1()->GetMean()); for (unsigned int i = 0; i < nParams; i++) fitResults.push_back(0.0); - goodData = MyFit(&fitResults); + if (internalDebug == true) + { + cout << "[Vx3DHLTAnalyzer]::\t@@@ Beam Spot parameters - prefit @@@" << endl; + + printFitParams(fitResults); + + cout << "Runnumber " << runNumber << endl; + cout << "BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) << " " << (beginTimeOfFit >> 32) << endl; + cout << "EndTimeOfFit " << formatTime(endTimeOfFit >> 32) << " " << (endTimeOfFit >> 32) << endl; + cout << "LumiRange " << beginLumiOfFit << " - " << endLumiOfFit << endl; + } + + goodData = MyFit(&fitResults); if (internalDebug == true) { + cout << "[Vx3DHLTAnalyzer]::\t@@@ Beam Spot parameters - postfit @@@" << endl; + + printFitParams(fitResults); + cout << "goodData --> " << goodData << endl; cout << "Used vertices --> " << counterVx << endl; - cout << "var x --> " << fitResults[0] << " +/- " << fitResults[0+nParams] << endl; - cout << "var y --> " << fitResults[1] << " +/- " << fitResults[1+nParams] << endl; - cout << "var z --> " << fitResults[2] << " +/- " << fitResults[2+nParams] << endl; - cout << "cov xy --> " << fitResults[3] << " +/- " << fitResults[3+nParams] << endl; - cout << "dydz --> " << fitResults[4] << " +/- " << fitResults[4+nParams] << endl; - cout << "dxdz --> " << fitResults[5] << " +/- " << fitResults[5+nParams] << endl; - cout << "mean x --> " << fitResults[6] << " +/- " << fitResults[6+nParams] << endl; - cout << "mean y --> " << fitResults[7] << " +/- " << fitResults[7+nParams] << endl; - cout << "mean z --> " << fitResults[8] << " +/- " << fitResults[8+nParams] << endl; } if (goodData == 0) @@ -856,7 +950,7 @@ void Vx3DHLTAnalyzer::endLuminosityBlock (const LuminosityBlock& lumiBlock, cons vals.push_back(std::pow(std::fabs(fitResults[0+nParams]) / (2.*std::sqrt(std::fabs(fitResults[0]))),2.)); vals.push_back(std::pow(std::fabs(fitResults[1+nParams]) / (2.*std::sqrt(std::fabs(fitResults[1]))),2.)); } - else for (unsigned int i = 0; i < 8*2; i++) vals.push_back(0.0); + else for (unsigned int i = 0; i < (nParams-1)*2; i++) vals.push_back(0.0); fitResults.clear(); } @@ -889,7 +983,7 @@ void Vx3DHLTAnalyzer::endLuminosityBlock (const LuminosityBlock& lumiBlock, cons else { goodData = -2; - for (unsigned int i = 0; i < 8*2; i++) vals.push_back(0.0); + for (unsigned int i = 0; i < (nParams-1)*2; i++) vals.push_back(0.0); } } @@ -911,49 +1005,39 @@ void Vx3DHLTAnalyzer::endLuminosityBlock (const LuminosityBlock& lumiBlock, cons // vals[14] = err^2 BeamWidthX // vals[15] = err^2 BeamWidthY - // "goodData" CODE: - // 0 == OK --> Reset - // -2 == NO OK - not enough "minNentries" --> Wait for more lumisections - // Any other number == NO OK --> Reset - numberFits++; + writeToFile(&vals, beginTimeOfFit, endTimeOfFit, beginLumiOfFit, endLumiOfFit, 3); + if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tUsed vertices: " << counterVx << endl; + + statusCounter->getTH1()->SetBinContent(lastLumiOfFit, (double)goodData); + statusCounter->getTH1()->SetBinError(lastLumiOfFit, 1e-3); + if (goodData == 0) { - writeToFile(&vals, beginTimeOfFit, endTimeOfFit, beginLumiOfFit, endLumiOfFit, 3); - if ((internalDebug == true) && (outputDebugFile.is_open() == true)) outputDebugFile << "Used vertices: " << counterVx << endl; - numberGoodFits++; - histTitle << "Fitted Beam Spot [cm] (Lumi start: " << beginLumiOfFit << " - Lumi end: " << endLumiOfFit << ")"; - if (lumiCounterHisto >= maxLumiIntegration) reset("whole"); - else reset("partial"); + histTitle << "Ongoing: fitted lumis " << beginLumiOfFit << " - " << endLumiOfFit; + reset("whole"); } else { - writeToFile(&vals, beginTimeOfFit, endTimeOfFit, beginLumiOfFit, endLumiOfFit, -1); - if ((internalDebug == true) && (outputDebugFile.is_open() == true)) outputDebugFile << "Used vertices: " << counterVx << endl; - - if (goodData == -2) - { - histTitle << "Fitted Beam Spot [cm] (not enough statistics)"; - if (lumiCounter >= maxLumiIntegration) reset("whole"); - else reset("hitCounter"); - } - else + if (goodData == -2) histTitle << "Ongoing: not enough evts (" << lumiCounter << " - " << maxLumiIntegration << " lumis)"; + else histTitle << "Ongoing: temporary problems (" << lumiCounter << " - " << maxLumiIntegration << " lumis)"; + + if (lumiCounter > maxLumiIntegration) { - histTitle << "Fitted Beam Spot [cm] (problems)"; - if (lumiCounterHisto >= maxLumiIntegration) reset("whole"); - else reset("partial"); - - counterVx = 0; + statusCounter->getTH1()->SetBinContent(lastLumiOfFit, -5); + statusCounter->getTH1()->SetBinError(lastLumiOfFit, 1e-3); + reset("whole"); } + else reset("hitCounter"); } - reportSummary->Fill(numberFits != 0 ? (double)numberGoodFits/(double)numberFits : 0.0); - reportSummaryMap->Fill(0.5, 0.5, numberFits != 0 ? (double)numberGoodFits/(double)numberFits : 0.0); + reportSummary->Fill((numberFits != 0 ? ((double)numberGoodFits) / ((double)numberFits) : -1)); + reportSummaryMap->getTH1()->SetBinContent(1, 1, (numberFits != 0 ? ((double)numberGoodFits) / ((double)numberFits) : -1)); fitResults->setAxisTitle(histTitle.str().c_str(), 1); - + fitResults->setBinContent(1, 9, vals[0]); fitResults->setBinContent(1, 8, vals[1]); fitResults->setBinContent(1, 7, vals[2]); @@ -978,229 +1062,218 @@ void Vx3DHLTAnalyzer::endLuminosityBlock (const LuminosityBlock& lumiBlock, cons TF1* myLinFit = new TF1("myLinFit", "[0] + [1]*x", mXlumi->getTH1()->GetXaxis()->GetXmin(), mXlumi->getTH1()->GetXaxis()->GetXmax()); myLinFit->SetLineColor(2); myLinFit->SetLineWidth(2); - myLinFit->SetParName(0,"Intercept"); + myLinFit->SetParName(0,"Inter."); myLinFit->SetParName(1,"Slope"); - mXlumi->ShiftFillLast(vals[0], std::sqrt(vals[8]), nLumiReset); + mXlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[0]); + mXlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[8])); myLinFit->SetParameter(0, mXlumi->getTH1()->GetMean(2)); myLinFit->SetParameter(1, 0.0); mXlumi->getTH1()->Fit(myLinFit,"QR"); - mYlumi->ShiftFillLast(vals[1], std::sqrt(vals[9]), nLumiReset); + mYlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[1]); + mYlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[9])); myLinFit->SetParameter(0, mYlumi->getTH1()->GetMean(2)); myLinFit->SetParameter(1, 0.0); mYlumi->getTH1()->Fit(myLinFit,"QR"); - mZlumi->ShiftFillLast(vals[2], std::sqrt(vals[10]), nLumiReset); + mZlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[2]); + mZlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[10])); myLinFit->SetParameter(0, mZlumi->getTH1()->GetMean(2)); myLinFit->SetParameter(1, 0.0); mZlumi->getTH1()->Fit(myLinFit,"QR"); - sXlumi->ShiftFillLast(vals[6], std::sqrt(vals[14]), nLumiReset); + sXlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[6]); + sXlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[14])); myLinFit->SetParameter(0, sXlumi->getTH1()->GetMean(2)); myLinFit->SetParameter(1, 0.0); sXlumi->getTH1()->Fit(myLinFit,"QR"); - sYlumi->ShiftFillLast(vals[7], std::sqrt(vals[15]), nLumiReset); + sYlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[7]); + sYlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[15])); myLinFit->SetParameter(0, sYlumi->getTH1()->GetMean(2)); myLinFit->SetParameter(1, 0.0); sYlumi->getTH1()->Fit(myLinFit,"QR"); - sZlumi->ShiftFillLast(vals[3], std::sqrt(vals[11]), nLumiReset); + sZlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[3]); + sZlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[11])); myLinFit->SetParameter(0, sZlumi->getTH1()->GetMean(2)); myLinFit->SetParameter(1, 0.0); sZlumi->getTH1()->Fit(myLinFit,"QR"); - dxdzlumi->ShiftFillLast(vals[4], std::sqrt(vals[12]), nLumiReset); + dxdzlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[4]); + dxdzlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[12])); myLinFit->SetParameter(0, dxdzlumi->getTH1()->GetMean(2)); myLinFit->SetParameter(1, 0.0); dxdzlumi->getTH1()->Fit(myLinFit,"QR"); - dydzlumi->ShiftFillLast(vals[5], std::sqrt(vals[13]), nLumiReset); + dydzlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[5]); + dydzlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[13])); myLinFit->SetParameter(0, dydzlumi->getTH1()->GetMean(2)); myLinFit->SetParameter(1, 0.0); dydzlumi->getTH1()->Fit(myLinFit,"QR"); - goodVxCounter->ShiftFillLast((double)counterVx, std::sqrt((double)counterVx), nLumiReset); - myLinFit->SetParameter(0, goodVxCounter->getTH1()->GetMean(2)); - myLinFit->SetParameter(1, 0.0); - goodVxCounter->getTH1()->Fit(myLinFit,"QR"); - delete myLinFit; // Exponential fit to the historical plot - TF1* myExpFit = new TF1("myExpFit", "[0]*exp(-x/[1])", hitCountHistory->getTH1()->GetXaxis()->GetXmin(), hitCountHistory->getTH1()->GetXaxis()->GetXmax()); + TF1* myExpFit = new TF1("myExpFit", "[0]*exp(-x/[1])", hitCounter->getTH1()->GetXaxis()->GetXmin(), hitCounter->getTH1()->GetXaxis()->GetXmax()); myExpFit->SetLineColor(2); myExpFit->SetLineWidth(2); - myExpFit->SetParName(0,"Amplitude"); + myExpFit->SetParName(0,"Ampli."); myExpFit->SetParName(1,"#tau"); - myExpFit->SetParameter(0, hitCountHistory->getTH1()->GetBinContent(1)); - myExpFit->SetParameter(1, nBinsWholeHistory/2); - hitCountHistory->getTH1()->Fit(myExpFit,"QR"); + myExpFit->SetParameter(0, hitCounter->getTH1()->GetMaximum()); + myExpFit->SetParameter(1, nLumiXaxisRange/2); + hitCounter->getTH1()->Fit(myExpFit,"QR"); - if (lastLumiOfFit % prescaleHistory == 0) - { - goodVxCountHistory->getTH1()->SetBinContent(lastLumiOfFit, (double)counterVx); - goodVxCountHistory->getTH1()->SetBinError(lastLumiOfFit, std::sqrt((double)counterVx)); - - myExpFit->SetParameter(0, goodVxCountHistory->getTH1()->GetBinContent(1)); - myExpFit->SetParameter(1, nBinsWholeHistory/2); - goodVxCountHistory->getTH1()->Fit(myExpFit,"QR"); - } + goodVxCounter->getTH1()->SetBinContent(lastLumiOfFit, (double)counterVx); + goodVxCounter->getTH1()->SetBinError(lastLumiOfFit, std::sqrt((double)counterVx)); + + myExpFit->SetParameter(0, goodVxCounter->getTH1()->GetMaximum()); + myExpFit->SetParameter(1, nLumiXaxisRange/2); + goodVxCounter->getTH1()->Fit(myExpFit,"QR"); delete myExpFit; - vals.clear(); } - else if (nLumiReset == 0) + else if ((nLumiFit != 0) && (lumiCounter%nLumiFit != 0) && (beginTimeOfFit != 0) && (runNumber != 0)) { - histTitle << "Fitted Beam Spot [cm] (no ongoing fits)"; + histTitle << "Ongoing: accumulating evts (" << lumiCounter%nLumiFit << " - " << nLumiFit << " in " << lumiCounter << " - " << maxLumiIntegration << " lumis)"; fitResults->setAxisTitle(histTitle.str().c_str(), 1); - hitCounter->ShiftFillLast(totalHits, std::sqrt(totalHits), 1); - reset("nohisto"); + if ((debugMode == true) && (outputDebugFile.is_open() == true)) + { + outputDebugFile << "Runnumber " << runNumber << endl; + outputDebugFile << "BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) << " " << (beginTimeOfFit >> 32) << endl; + outputDebugFile << "BeginLumiRange " << beginLumiOfFit << endl; + outputDebugFile << histTitle.str().c_str() << "\n" << endl; + } } -} + else if ((nLumiFit == 0) || (beginTimeOfFit == 0) || (runNumber == 0)) + { + histTitle << "Ongoing: no ongoing fits"; + fitResults->setAxisTitle(histTitle.str().c_str(), 1); + if ((debugMode == true) && (outputDebugFile.is_open() == true)) outputDebugFile << histTitle.str().c_str() << "\n" << endl; + endLumiOfFit = lumiBlock.luminosityBlock(); -void Vx3DHLTAnalyzer::beginJob () -{ - // ### Set internal variables ### - prescaleHistory = 1; // Set the number of lumis to update historical plot - maxLumiIntegration = 15; // If failing fits, this is the maximum number of integrated lumis after which a reset is issued - minVxDoF = 10.; // Good vertex selection cut - // For vertex fitter without track-weight: d.o.f. = 2*NTracks - 3 - // For vertex fitter with track-weight: d.o.f. = sum_NTracks(2*track_weight) - 3 - internalDebug = false; - considerVxCovariance = true; // Deconvolute vertex covariance matrix - pi = 3.141592653589793238; - // ############################## -} + hitCounter->getTH1()->SetBinContent(endLumiOfFit, (double)totalHits); + hitCounter->getTH1()->SetBinError(endLumiOfFit, std::sqrt((double)totalHits)); + reset("whole"); + } -void Vx3DHLTAnalyzer::endJob () { reset("scratch"); } + if (internalDebug == true) cout << "[Vx3DHLTAnalyzer]::\tHistogram title: " << histTitle.str() << endl; +} -void Vx3DHLTAnalyzer::beginRun (const Run& iRun, const EventSetup& iSetup) +void Vx3DHLTAnalyzer::bookHistograms(DQMStore::IBooker & ibooker, Run const & iRun, EventSetup const & /* iSetup */) { - DQMStore* dbe = 0; - dbe = Service().operator->(); - - // ### Set internal variables ### - nBinsHistoricalPlot = 80; - nBinsWholeHistory = 3000; // Correspond to about 20h of data taking: 20h * 60min * 60s / 23s per lumi-block = 3130 - // ############################## - - if ( dbe ) - { - dbe->setCurrentFolder("BeamPixel"); - - Vx_X = dbe->book1D("vertex x", "Primary Vertex X Coordinate Distribution", int(rint(xRange/xStep)), -xRange/2., xRange/2.); - Vx_Y = dbe->book1D("vertex y", "Primary Vertex Y Coordinate Distribution", int(rint(yRange/yStep)), -yRange/2., yRange/2.); - Vx_Z = dbe->book1D("vertex z", "Primary Vertex Z Coordinate Distribution", int(rint(zRange/zStep)), -zRange/2., zRange/2.); - Vx_X->setAxisTitle("Primary Vertices X [cm]",1); - Vx_X->setAxisTitle("Entries [#]",2); - Vx_Y->setAxisTitle("Primary Vertices Y [cm]",1); - Vx_Y->setAxisTitle("Entries [#]",2); - Vx_Z->setAxisTitle("Primary Vertices Z [cm]",1); - Vx_Z->setAxisTitle("Entries [#]",2); - - mXlumi = dbe->book1D("muX vs lumi", "#mu_{x} vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5); - mYlumi = dbe->book1D("muY vs lumi", "#mu_{y} vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5); - mZlumi = dbe->book1D("muZ vs lumi", "#mu_{z} vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5); - mXlumi->setAxisTitle("Lumisection [#]",1); - mXlumi->setAxisTitle("#mu_{x} [cm]",2); - mXlumi->getTH1()->SetOption("E1"); - mYlumi->setAxisTitle("Lumisection [#]",1); - mYlumi->setAxisTitle("#mu_{y} [cm]",2); - mYlumi->getTH1()->SetOption("E1"); - mZlumi->setAxisTitle("Lumisection [#]",1); - mZlumi->setAxisTitle("#mu_{z} [cm]",2); - mZlumi->getTH1()->SetOption("E1"); - - sXlumi = dbe->book1D("sigmaX vs lumi", "#sigma_{x} vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5); - sYlumi = dbe->book1D("sigmaY vs lumi", "#sigma_{y} vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5); - sZlumi = dbe->book1D("sigmaZ vs lumi", "#sigma_{z} vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5); - sXlumi->setAxisTitle("Lumisection [#]",1); - sXlumi->setAxisTitle("#sigma_{x} [cm]",2); - sXlumi->getTH1()->SetOption("E1"); - sYlumi->setAxisTitle("Lumisection [#]",1); - sYlumi->setAxisTitle("#sigma_{y} [cm]",2); - sYlumi->getTH1()->SetOption("E1"); - sZlumi->setAxisTitle("Lumisection [#]",1); - sZlumi->setAxisTitle("#sigma_{z} [cm]",2); - sZlumi->getTH1()->SetOption("E1"); - - dxdzlumi = dbe->book1D("dxdz vs lumi", "dX/dZ vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5); - dydzlumi = dbe->book1D("dydz vs lumi", "dY/dZ vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5); - dxdzlumi->setAxisTitle("Lumisection [#]",1); - dxdzlumi->setAxisTitle("dX/dZ [rad]",2); - dxdzlumi->getTH1()->SetOption("E1"); - dydzlumi->setAxisTitle("Lumisection [#]",1); - dydzlumi->setAxisTitle("dY/dZ [rad]",2); - dydzlumi->getTH1()->SetOption("E1"); - - Vx_ZX = dbe->book2D("vertex zx", "Primary Vertex ZX Coordinate Distribution", int(rint(zRange/zStep/5.)), -zRange/2., zRange/2., int(rint(xRange/xStep/5.)), -xRange/2., xRange/2.); - Vx_ZY = dbe->book2D("vertex zy", "Primary Vertex ZY Coordinate Distribution", int(rint(zRange/zStep/5.)), -zRange/2., zRange/2., int(rint(yRange/yStep/5.)), -yRange/2., yRange/2.); - Vx_XY = dbe->book2D("vertex xy", "Primary Vertex XY Coordinate Distribution", int(rint(xRange/xStep/5.)), -xRange/2., xRange/2., int(rint(yRange/yStep/5.)), -yRange/2., yRange/2.); - Vx_ZX->setAxisTitle("Primary Vertices Z [cm]",1); - Vx_ZX->setAxisTitle("Primary Vertices X [cm]",2); - Vx_ZX->setAxisTitle("Entries [#]",3); - Vx_ZY->setAxisTitle("Primary Vertices Z [cm]",1); - Vx_ZY->setAxisTitle("Primary Vertices Y [cm]",2); - Vx_ZY->setAxisTitle("Entries [#]",3); - Vx_XY->setAxisTitle("Primary Vertices X [cm]",1); - Vx_XY->setAxisTitle("Primary Vertices Y [cm]",2); - Vx_XY->setAxisTitle("Entries [#]",3); - - hitCounter = dbe->book1D("pixelHits vs lumi", "# Pixel-Hits vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5); - hitCounter->setAxisTitle("Lumisection [#]",1); - hitCounter->setAxisTitle("Pixel-Hits [#]",2); - hitCounter->getTH1()->SetOption("E1"); - - hitCountHistory = dbe->book1D("hist pixelHits vs lumi", "History: # Pixel-Hits vs. Lumi", nBinsWholeHistory, 0.5, (double)nBinsWholeHistory+0.5); - hitCountHistory->setAxisTitle("Lumisection [#]",1); - hitCountHistory->setAxisTitle("Pixel-Hits [#]",2); - hitCountHistory->getTH1()->SetOption("E1"); - - goodVxCounter = dbe->book1D("good vertices vs lumi", "# Good vertices vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5); - goodVxCounter->setAxisTitle("Lumisection [#]",1); - goodVxCounter->setAxisTitle("Good vertices [#]",2); - goodVxCounter->getTH1()->SetOption("E1"); - - goodVxCountHistory = dbe->book1D("hist good vx vs lumi", "History: # Good vx vs. Lumi", nBinsWholeHistory, 0.5, (double)nBinsWholeHistory+0.5); - goodVxCountHistory->setAxisTitle("Lumisection [#]",1); - goodVxCountHistory->setAxisTitle("Good vertices [#]",2); - goodVxCountHistory->getTH1()->SetOption("E1"); - - fitResults = dbe->book2D("fit results","Results of Beam Spot Fit", 2, 0., 2., 9, 0., 9.); - fitResults->setAxisTitle("Fitted Beam Spot [cm]", 1); - fitResults->setBinLabel(9, "X", 2); - fitResults->setBinLabel(8, "Y", 2); - fitResults->setBinLabel(7, "Z", 2); - fitResults->setBinLabel(6, "#sigma_{Z}", 2); - fitResults->setBinLabel(5, "#frac{dX}{dZ}[rad]", 2); - fitResults->setBinLabel(4, "#frac{dY}{dZ}[rad]", 2); - fitResults->setBinLabel(3, "#sigma_{X}", 2); - fitResults->setBinLabel(2, "#sigma_{Y}", 2); - fitResults->setBinLabel(1, "Vertices", 2); - fitResults->setBinLabel(1, "Value", 1); - fitResults->setBinLabel(2, "Stat. Error", 1); - fitResults->getTH1()->SetOption("text"); - - dbe->setCurrentFolder("BeamPixel/EventInfo"); - reportSummary = dbe->bookFloat("reportSummary"); - reportSummary->Fill(-1.); - reportSummaryMap = dbe->book2D("reportSummaryMap","Pixel-Vertices Beam Spot: % Good Fits", 1, 0., 1., 1, 0., 1.); - reportSummaryMap->Fill(0.5, 0.5, -1.); - dbe->setCurrentFolder("BeamPixel/EventInfo/reportSummaryContents"); - - // Convention for reportSummary and reportSummaryMap: - // - 0% at the moment of creation of the histogram - // - n% numberGoodFits / numberFits + ibooker.setCurrentFolder("BeamPixel"); + + Vx_X = ibooker.book1D("F - vertex x", "Primary Vertex X Coordinate Distribution", int(rint(xRange/xStep)), -xRange/2., xRange/2.); + Vx_Y = ibooker.book1D("F - vertex y", "Primary Vertex Y Coordinate Distribution", int(rint(yRange/yStep)), -yRange/2., yRange/2.); + Vx_Z = ibooker.book1D("F - vertex z", "Primary Vertex Z Coordinate Distribution", int(rint(zRange/zStep)), -zRange/2., zRange/2.); + Vx_X->setAxisTitle("Primary Vertices X [cm]",1); + Vx_X->setAxisTitle("Entries [#]",2); + Vx_Y->setAxisTitle("Primary Vertices Y [cm]",1); + Vx_Y->setAxisTitle("Entries [#]",2); + Vx_Z->setAxisTitle("Primary Vertices Z [cm]",1); + Vx_Z->setAxisTitle("Entries [#]",2); + + mXlumi = ibooker.book1D("B - muX vs lumi", "#mu_{x} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange)+0.5); + mYlumi = ibooker.book1D("B - muY vs lumi", "#mu_{y} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange)+0.5); + mZlumi = ibooker.book1D("B - muZ vs lumi", "#mu_{z} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange)+0.5); + mXlumi->setAxisTitle("Lumisection [#]",1); + mXlumi->setAxisTitle("#mu_{x} [cm]",2); + mXlumi->getTH1()->SetOption("E1"); + mYlumi->setAxisTitle("Lumisection [#]",1); + mYlumi->setAxisTitle("#mu_{y} [cm]",2); + mYlumi->getTH1()->SetOption("E1"); + mZlumi->setAxisTitle("Lumisection [#]",1); + mZlumi->setAxisTitle("#mu_{z} [cm]",2); + mZlumi->getTH1()->SetOption("E1"); + + sXlumi = ibooker.book1D("C - sigmaX vs lumi", "#sigma_{x} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange)+0.5); + sYlumi = ibooker.book1D("C - sigmaY vs lumi", "#sigma_{y} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange)+0.5); + sZlumi = ibooker.book1D("C - sigmaZ vs lumi", "#sigma_{z} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange)+0.5); + sXlumi->setAxisTitle("Lumisection [#]",1); + sXlumi->setAxisTitle("#sigma_{x} [cm]",2); + sXlumi->getTH1()->SetOption("E1"); + sYlumi->setAxisTitle("Lumisection [#]",1); + sYlumi->setAxisTitle("#sigma_{y} [cm]",2); + sYlumi->getTH1()->SetOption("E1"); + sZlumi->setAxisTitle("Lumisection [#]",1); + sZlumi->setAxisTitle("#sigma_{z} [cm]",2); + sZlumi->getTH1()->SetOption("E1"); + + dxdzlumi = ibooker.book1D("D - dxdz vs lumi", "dX/dZ vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange)+0.5); + dydzlumi = ibooker.book1D("D - dydz vs lumi", "dY/dZ vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange)+0.5); + dxdzlumi->setAxisTitle("Lumisection [#]",1); + dxdzlumi->setAxisTitle("dX/dZ [rad]",2); + dxdzlumi->getTH1()->SetOption("E1"); + dydzlumi->setAxisTitle("Lumisection [#]",1); + dydzlumi->setAxisTitle("dY/dZ [rad]",2); + dydzlumi->getTH1()->SetOption("E1"); + + Vx_ZX = ibooker.book2D("E - vertex zx", "Primary Vertex ZX Coordinate Distribution", int(rint(zRange/zStep)), -zRange/2., zRange/2., int(rint(xRange/xStep)), -xRange/2., xRange/2.); + Vx_ZY = ibooker.book2D("E - vertex zy", "Primary Vertex ZY Coordinate Distribution", int(rint(zRange/zStep)), -zRange/2., zRange/2., int(rint(yRange/yStep)), -yRange/2., yRange/2.); + Vx_XY = ibooker.book2D("E - vertex xy", "Primary Vertex XY Coordinate Distribution", int(rint(xRange/xStep)), -xRange/2., xRange/2., int(rint(yRange/yStep)), -yRange/2., yRange/2.); + Vx_ZX->setAxisTitle("Primary Vertices Z [cm]",1); + Vx_ZX->setAxisTitle("Primary Vertices X [cm]",2); + Vx_ZX->setAxisTitle("Entries [#]",3); + Vx_ZY->setAxisTitle("Primary Vertices Z [cm]",1); + Vx_ZY->setAxisTitle("Primary Vertices Y [cm]",2); + Vx_ZY->setAxisTitle("Entries [#]",3); + Vx_XY->setAxisTitle("Primary Vertices X [cm]",1); + Vx_XY->setAxisTitle("Primary Vertices Y [cm]",2); + Vx_XY->setAxisTitle("Entries [#]",3); + + hitCounter = ibooker.book1D("H - pixelHits vs lumi", "# Pixel-Hits vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange)+0.5); + hitCounter->setAxisTitle("Lumisection [#]",1); + hitCounter->setAxisTitle("Pixel-Hits [#]",2); + hitCounter->getTH1()->SetOption("E1"); + + goodVxCounter = ibooker.book1D("G - good vertices vs lumi", "# Good vertices vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange)+0.5); + goodVxCounter->setAxisTitle("Lumisection [#]",1); + goodVxCounter->setAxisTitle("Good vertices [#]",2); + goodVxCounter->getTH1()->SetOption("E1"); - reset("scratch"); // Initialize histograms after creation - } + statusCounter = ibooker.book1D("I - app status vs lumi", "Status vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange)+0.5); + statusCounter->setAxisTitle("Lumisection [#]",1); + statusCounter->setAxisTitle("App. status [0 = OK]",2); + statusCounter->getTH1()->SetOption("E1"); + + fitResults = ibooker.book2D("A - fit results","Results of Beam Spot Fit", 2, 0., 2., 9, 0., 9.); + fitResults->setAxisTitle("Ongoing: bootstrapping", 1); + fitResults->setBinLabel(9, "X[cm]", 2); + fitResults->setBinLabel(8, "Y[cm]", 2); + fitResults->setBinLabel(7, "Z[cm]", 2); + fitResults->setBinLabel(6, "#sigma_{Z}[cm]", 2); + fitResults->setBinLabel(5, "#frac{dX}{dZ}[rad]", 2); + fitResults->setBinLabel(4, "#frac{dY}{dZ}[rad]", 2); + fitResults->setBinLabel(3, "#sigma_{X}[cm]", 2); + fitResults->setBinLabel(2, "#sigma_{Y}[cm]", 2); + fitResults->setBinLabel(1, "Vtx[#]", 2); + fitResults->setBinLabel(1, "Value", 1); + fitResults->setBinLabel(2, "Error (stat)", 1); + fitResults->getTH1()->SetOption("text"); + + + ibooker.setCurrentFolder("BeamPixel/EventInfo"); + + reportSummary = ibooker.bookFloat("reportSummary"); + reportSummary->Fill(-1); + reportSummaryMap = ibooker.book2D("reportSummaryMap","Pixel-Vertices Beam Spot: % Good Fits", 1, 0., 1., 1, 0., 1.); + reportSummaryMap->getTH1()->SetBinContent(1, 1, -1); + + ibooker.setCurrentFolder("BeamPixel/EventInfo/reportSummaryContents"); + + // Convention for reportSummary and reportSummaryMap: + // - -1% at the moment of creation of the histogram (i.e. white histogram) + // - n% numberGoodFits / numberFits + + + reset("scratch"); // Initialize histograms after creation } diff --git a/DQM/BeamMonitor/plugins/Vx3DHLTAnalyzer.h b/DQM/BeamMonitor/plugins/Vx3DHLTAnalyzer.h index 77746b4bcb110..4f763ce6dae6d 100644 --- a/DQM/BeamMonitor/plugins/Vx3DHLTAnalyzer.h +++ b/DQM/BeamMonitor/plugins/Vx3DHLTAnalyzer.h @@ -2,24 +2,18 @@ #define Vx3DHLTAnalyzer_H // -*- C++ -*- -// -// Package: Vx3DHLTAnalyzer -// Class: Vx3DHLTAnalyzer -// -/**\class Vx3DHLTAnalyzer Vx3DHLTAnalyzer.cc interface/Vx3DHLTAnalyzer.h - - Description: beam-spot monitor entirely based on pixel detector information - Implementation: the monitoring is based on a 3D fit to the vertex cloud +// Package: Vx3DHLTAnalyzer +// Class: Vx3DHLTAnalyzer + +/* + Class Vx3DHLTAnalyzer Vx3DHLTAnalyzer.cc plugins/Vx3DHLTAnalyzer.cc + + Description: beam-spot monitor entirely based on pixel detector information + Implementation: the monitoring is based on a 3D fit to the vertex cloud */ -// -// Original Author: Mauro Dinardo, 28 S-012, +41-22-767-8302, -// Created: Tue Feb 23 13:15:31 CET 2010 -// -*- C++ -*- -// -// Package: Vx3DHLTAnalyzer -// Class: Vx3DHLTAnalyzer -// +// Original Author: Mauro Dinardo, 28 S-012, +41-22-767-8302 +// Created: Tue Feb 23 13:15:31 CET 2010 #include @@ -30,6 +24,7 @@ #include "DQMServices/Core/interface/DQMStore.h" #include "DQMServices/Core/interface/MonitorElement.h" +#include "DQMServices/Core/interface/DQMEDAnalyzer.h" #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h" #include "DataFormats/VertexReco/interface/VertexFwd.h" @@ -42,11 +37,6 @@ #include -using namespace std; -using namespace reco; -using namespace edm; - - // ################# // # Fit variables # // ################# @@ -59,7 +49,7 @@ typedef struct double z; double Covariance[DIM][DIM]; } VertexType; -vector Vertices; +std::vector Vertices; bool considerVxCovariance; unsigned int counterVx; // Counts the number of vertices taken into account for the fit double maxTransRadius; // Max transverse radius in which the vertices must be [cm] @@ -72,104 +62,105 @@ double pi; double VxErrCorr; // Coefficient to compensate the under-estimation of the vertex errors -class Vx3DHLTAnalyzer : public EDAnalyzer { - public: - explicit Vx3DHLTAnalyzer (const ParameterSet&); - ~Vx3DHLTAnalyzer(); - - - private: - virtual void beginJob (); - virtual void analyze (const Event&, const EventSetup&); - virtual unsigned int HitCounter (const Event& iEvent); - virtual string formatTime (const time_t& t); - virtual int MyFit (vector* vals); - virtual void reset (string ResetType); - virtual void writeToFile (vector* vals, - TimeValue_t BeginTimeOfFit, - TimeValue_t EndTimeOfFit, - unsigned int BeginLumiOfFit, - unsigned int EndLumiOfFit, - int dataType); - virtual void beginLuminosityBlock (const LuminosityBlock& lumiBlock, const EventSetup& iSetup); - virtual void endLuminosityBlock (const LuminosityBlock& lumiBlock, const EventSetup& iSetup); - virtual void endJob (); - virtual void beginRun (const Run& iRun, const EventSetup& iSetup); - - - // ####################### - // # cfg file parameters # - // ####################### - EDGetTokenT vertexCollection; - EDGetTokenT pixelHitCollection; - bool debugMode; - unsigned int nLumiReset; - bool dataFromFit; - unsigned int minNentries; - double xRange; - double xStep; - double yRange; - double yStep; - double zRange; - double zStep; - string fileName; - - - // ############## - // # Histograms # - // ############## - MonitorElement* mXlumi; - MonitorElement* mYlumi; - MonitorElement* mZlumi; - - MonitorElement* sXlumi; - MonitorElement* sYlumi; - MonitorElement* sZlumi; - - MonitorElement* dxdzlumi; - MonitorElement* dydzlumi; - - MonitorElement* Vx_X; - MonitorElement* Vx_Y; - MonitorElement* Vx_Z; - - MonitorElement* Vx_ZX; - MonitorElement* Vx_ZY; - MonitorElement* Vx_XY; - - MonitorElement* goodVxCounter; - MonitorElement* goodVxCountHistory; - - MonitorElement* hitCounter; - MonitorElement* hitCountHistory; - - MonitorElement* reportSummary; - MonitorElement* reportSummaryMap; - - MonitorElement* fitResults; - - // ###################### - // # Internal variables # - // ###################### - ofstream outputFile; - ofstream outputDebugFile; - TimeValue_t beginTimeOfFit; - TimeValue_t endTimeOfFit; - unsigned int nBinsHistoricalPlot; - unsigned int nBinsWholeHistory; - unsigned int runNumber; - unsigned int lumiCounter; - unsigned int lumiCounterHisto; - unsigned int totalHits; - unsigned int maxLumiIntegration; - unsigned int prescaleHistory; - unsigned int numberGoodFits; - unsigned int numberFits; - unsigned int beginLumiOfFit; - unsigned int endLumiOfFit; - unsigned int lastLumiOfFit; - double minVxDoF; - bool internalDebug; +class Vx3DHLTAnalyzer : public DQMEDAnalyzer +{ + + + public: + Vx3DHLTAnalyzer (const edm::ParameterSet&); + ~Vx3DHLTAnalyzer (); + + + private: + void analyze (const edm::Event& iEvent, const edm::EventSetup& iSetup); + void beginLuminosityBlock (const edm::LuminosityBlock& lumiBlock, const edm::EventSetup& iSetup); + void endLuminosityBlock (const edm::LuminosityBlock& lumiBlock, const edm::EventSetup& iSetup); + void bookHistograms (DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override; + + unsigned int HitCounter (const edm::Event& iEvent); + std::string formatTime (const time_t& t); + int MyFit (std::vector* vals); + void reset (std::string ResetType); + void writeToFile (std::vector* vals, + edm::TimeValue_t BeginTimeOfFit, + edm::TimeValue_t EndTimeOfFit, + unsigned int BeginLumiOfFit, + unsigned int EndLumiOfFit, + int dataType); + void printFitParams (const std::vector& fitResults); + + + // ####################### + // # cfg file parameters # + // ####################### + edm::EDGetTokenT vertexCollection; + edm::EDGetTokenT pixelHitCollection; + bool debugMode; + bool dataFromFit; + unsigned int nLumiFit; + unsigned int maxLumiIntegration; + unsigned int nLumiXaxisRange; + unsigned int minNentries; + double xRange; + double xStep; + double yRange; + double yStep; + double zRange; + double zStep; + double minVxDoF; + double minVxWgt; + std::string fileName; + + + // ############## + // # Histograms # + // ############## + MonitorElement* mXlumi; + MonitorElement* mYlumi; + MonitorElement* mZlumi; + + MonitorElement* sXlumi; + MonitorElement* sYlumi; + MonitorElement* sZlumi; + + MonitorElement* dxdzlumi; + MonitorElement* dydzlumi; + + MonitorElement* Vx_X; + MonitorElement* Vx_Y; + MonitorElement* Vx_Z; + + MonitorElement* Vx_ZX; + MonitorElement* Vx_ZY; + MonitorElement* Vx_XY; + + MonitorElement* goodVxCounter; + MonitorElement* hitCounter; + MonitorElement* statusCounter; + + MonitorElement* reportSummary; + MonitorElement* reportSummaryMap; + + MonitorElement* fitResults; + + + // ###################### + // # Internal variables # + // ###################### + std::ofstream outputFile; + std::ofstream outputDebugFile; + edm::TimeValue_t beginTimeOfFit; + edm::TimeValue_t endTimeOfFit; + unsigned int runNumber; + unsigned int lumiCounter; + unsigned int totalHits; + unsigned int numberGoodFits; + unsigned int numberFits; + unsigned int beginLumiOfFit; + unsigned int endLumiOfFit; + unsigned int lastLumiOfFit; + unsigned int nParams; + bool internalDebug; }; #endif diff --git a/DQM/Integration/python/clients/beampixel_dqm_sourceclient-live_cfg.py b/DQM/Integration/python/clients/beampixel_dqm_sourceclient-live_cfg.py index 56578a4fe1f27..ed5e92673c07d 100644 --- a/DQM/Integration/python/clients/beampixel_dqm_sourceclient-live_cfg.py +++ b/DQM/Integration/python/clients/beampixel_dqm_sourceclient-live_cfg.py @@ -6,14 +6,17 @@ #---------------------------- # Common for PP and HI running #---------------------------- -### @@@@@@ Comment when running locally @@@@@@ ### +# Use this to run locally (for testing purposes) +#process.load("DQM.Integration.config.fileinputsource_cfi") +# Otherwise use this process.load("DQM.Integration.config.inputsource_cfi") + #---------------------------- # HLT Filter #---------------------------- # 0=random, 1=physics, 2=calibration, 3=technical -process.hltTriggerTypeFilter = cms.EDFilter("HLTTriggerTypeFilter",SelectedTriggerType = cms.int32(1)) +process.hltTriggerTypeFilter = cms.EDFilter("HLTTriggerTypeFilter", SelectedTriggerType = cms.int32(1)) #---------------------------- @@ -30,14 +33,19 @@ process.load("Configuration.StandardSequences.GeometryRecoDB_cff") process.load('Configuration.StandardSequences.MagneticField_AutoFromDBCurrent_cff') process.load("Configuration.StandardSequences.RawToDigi_Data_cff") +# Use this to run locally (for testing purposes) +#process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_condDBv2_cff') +#from Configuration.AlCa.GlobalTag_condDBv2 import GlobalTag +#process.GlobalTag = GlobalTag(process.GlobalTag, '74X_dataRun2_Prompt_v0', '') +# Otherwise use this process.load("DQM.Integration.config.FrontierCondition_GT_cfi") #---------------------------- # Define Sequences #---------------------------- -process.dqmmodules = cms.Sequence(process.dqmEnv + process.dqmSaver) -process.phystrigger = cms.Sequence(process.hltTriggerTypeFilter) +process.dqmModules = cms.Sequence(process.dqmEnv + process.dqmSaver) +process.physTrigger = cms.Sequence(process.hltTriggerTypeFilter) #---------------------------- @@ -71,68 +79,69 @@ process.siPixelDigis.InputLabel = cms.InputTag("rawDataCollector") process.siStripDigis.ProductLabel = cms.InputTag("rawDataCollector") - process.load("Configuration.StandardSequences.Reconstruction_cff") + process.load('Configuration.StandardSequences.Reconstruction_Data_cff') #---------------------------- - # pixelVertexDQM Configuration + # pixelVertexDQM Config #---------------------------- process.pixelVertexDQM = cms.EDAnalyzer("Vx3DHLTAnalyzer", vertexCollection = cms.untracked.InputTag("pixelVertices"), - #pixelHitCollection = cms.untracked.InputTag("siPixelRecHits"), pixelHitCollection = cms.untracked.InputTag("siPixelRecHitsPreSplitting"), debugMode = cms.bool(True), - nLumiReset = cms.uint32(2), + nLumiFit = cms.uint32(2), + maxLumiIntegration = cms.uint32(15), + nLumiXaxisRange = cms.uint32(3000), dataFromFit = cms.bool(True), minNentries = cms.uint32(20), # If the histogram has at least "minNentries" then extract Mean and RMS, # or, if we are performing the fit, the number of vertices must be greater - # than minNentries otherwise it waits for other nLumiReset - xRange = cms.double(1.0), + # than minNentries otherwise it waits for other nLumiFit + xRange = cms.double(0.8), xStep = cms.double(0.001), - yRange = cms.double(1.0), + yRange = cms.double(0.8), yStep = cms.double(0.001), zRange = cms.double(30.0), - zStep = cms.double(0.05), - VxErrCorr = cms.double(1.3), # Keep checking this with later release - fileName = cms.string("/nfshome0/yumiceva/BeamMonitorDQM/BeamPixelResults.txt")) - if process.dqmRunConfig.type.value() is "playback": - process.pixelVertexDQM.fileName = cms.string("/nfshome0/dqmdev/BeamMonitorDQM/BeamPixelResults.txt") + zStep = cms.double(0.04), + VxErrCorr = cms.double(1.3), # Keep checking this with later release + minVxDoF = cms.double(10.0), + minVxWgt = cms.double(0.5), + fileName = cms.string("/nfshome0/dqmdev/BeamMonitorDQM/BeamPixelResults.txt")) + if process.dqmSaver.producer.value() is "Playback": + process.pixelVertexDQM.fileName = cms.string("/nfshome0/dqmdev/BeamMonitorDQM/BeamPixelResults.txt") else: - process.pixelVertexDQM.fileName = cms.string("/nfshome0/dqmpro/BeamMonitorDQM/BeamPixelResults.txt") - + process.pixelVertexDQM.fileName = cms.string("/nfshome0/dqmpro/BeamMonitorDQM/BeamPixelResults.txt") + print "[beampixel_dqm_sourceclient-live_cfg]::saving DIP file into " + str(process.pixelVertexDQM.fileName) - - process.load("RecoVertex.PrimaryVertexProducer.OfflinePixel3DPrimaryVertices_cfi") - #pixel track/vertices reco + #---------------------------- + # Pixel-Tracks&Vertices Config + #---------------------------- process.load("RecoPixelVertexing.Configuration.RecoPixelVertexing_cff") process.pixelVertices.TkFilterParameters.minPt = process.pixelTracks.RegionFactoryPSet.RegionPSet.ptMin process.offlinePrimaryVertices.TrackLabel = cms.InputTag("pixelTracks") - #process.dqmBeamMonitor.PVFitter.errorScale = 1.25 #keep checking this with new release expected close to 1.2 - from RecoTracker.TkSeedingLayers.PixelLayerTriplets_cfi import * process.PixelLayerTriplets.BPix.HitProducer = cms.string('siPixelRecHitsPreSplitting') process.PixelLayerTriplets.FPix.HitProducer = cms.string('siPixelRecHitsPreSplitting') from RecoPixelVertexing.PixelTrackFitting.PixelTracks_cff import * process.pixelTracks.OrderedHitsFactoryPSet.GeneratorPSet.SeedComparitorPSet.clusterShapeCacheSrc = cms.InputTag('siPixelClusterShapeCachePreSplitting') + #---------------------------- - # Pixel-Vertices Configuration + # Pixel-Tracks&Vertices Reco #---------------------------- - process.reconstruction_step = cms.Sequence(process.siPixelDigis* + process.reconstructionStep = cms.Sequence(process.siPixelDigis* process.offlineBeamSpot* process.siPixelClustersPreSplitting* process.siPixelRecHitsPreSplitting* process.siPixelClusterShapeCachePreSplitting* - process.recopixelvertexing - ) + process.recopixelvertexing) #---------------------------- # Define Path #---------------------------- - process.p = cms.Path(process.phystrigger*process.reconstruction_step*process.pixelVertexDQM*process.dqmmodules) + process.p = cms.Path(process.scalersRawToDigi*process.physTrigger*process.reconstructionStep*process.pixelVertexDQM*process.dqmModules) @@ -163,48 +172,49 @@ #---------------------------- - # pixelVertexDQM Configuration + # pixelVertexDQM Config #---------------------------- process.pixelVertexDQM = cms.EDAnalyzer("Vx3DHLTAnalyzer", vertexCollection = cms.untracked.InputTag("hiSelectedVertex"), pixelHitCollection = cms.untracked.InputTag("siPixelRecHits"), debugMode = cms.bool(True), - nLumiReset = cms.uint32(5), + nLumiFit = cms.uint32(5), + maxLumiIntegration = cms.uint32(15), + nLumiXaxisRange = cms.uint32(3000), dataFromFit = cms.bool(True), minNentries = cms.uint32(20), # If the histogram has at least "minNentries" then extract Mean and RMS, # or, if we are performing the fit, the number of vertices must be greater - # than minNentries otherwise it waits for other nLumiReset - xRange = cms.double(1.0), + # than minNentries otherwise it waits for other nLumiFit + xRange = cms.double(0.8), xStep = cms.double(0.001), - yRange = cms.double(1.0), + yRange = cms.double(0.8), yStep = cms.double(0.001), zRange = cms.double(30.0), - zStep = cms.double(0.05), - VxErrCorr = cms.double(1.3), # Keep checking this with later release - fileName = cms.string("/nfshome0/yumiceva/BeamMonitorDQM/BeamPixelResults.txt")) - if process.dqmRunConfig.type.value() is "playback": - process.pixelVertexDQM.fileName = cms.string("/nfshome0/dqmdev/BeamMonitorDQM/BeamPixelResults.txt") + zStep = cms.double(0.04), + VxErrCorr = cms.double(1.3), # Keep checking this with later release + minVxDoF = cms.double(10.0), + minVxWgt = cms.double(0.5), + fileName = cms.string("/nfshome0/dqmdev/BeamMonitorDQM/BeamPixelResults.txt")) + if process.dqmSaver.producer.value() is "Playback": + process.pixelVertexDQM.fileName = cms.string("/nfshome0/dqmdev/BeamMonitorDQM/BeamPixelResults.txt") else: - process.pixelVertexDQM.fileName = cms.string("/nfshome0/dqmpro/BeamMonitorDQM/BeamPixelResults.txt") + process.pixelVertexDQM.fileName = cms.string("/nfshome0/dqmpro/BeamMonitorDQM/BeamPixelResults.txt") + print "[beampixel_dqm_sourceclient-live_cfg]::saving DIP file into " + str(process.pixelVertexDQM.fileName) #---------------------------- - # Pixel-Vertices Configuration + # Pixel-Tracks&Vertices Reco #---------------------------- - process.reconstruction_step = cms.Sequence(process.siPixelDigis* - process.offlineBeamSpot* - process.siPixelClusters* - process.siPixelRecHits* - process.hiPixelVertices* - process.siPixelClusterShapeCache* - process.PixelLayerTriplets* - process.pixelTracks* - process.pixelVertices) + process.reconstructionStep = cms.Sequence(process.siPixelDigis* + process.offlineBeamSpot* + process.siPixelClusters* + process.siPixelRecHits* + process.hiPixelVertices* + process.hiPixel3PrimTracks) #---------------------------- # Define Path #---------------------------- - process.p = cms.Path(process.phystrigger*process.reconstruction_step*process.pixelVertexDQM*process.dqmmodules) - + process.p = cms.Path(process.scalersRawToDigi*process.physTrigger*process.reconstructionStep*process.pixelVertexDQM*process.dqmModules)