Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-04-29 23:10:16

0001 /*
0002   \File Vx3DHLTAnalyzer.cc
0003   \Display Beam-spot monitor entirely based on pixel detector information
0004            the monitoring is based on a 3D fit to the vertex cloud
0005   \Author Mauro Dinardo
0006   \Version $ Revision: 3.5 $
0007   \Date $ Date: 2010/23/02 13:15:00 $
0008 */
0009 
0010 #include "DQM/BeamMonitor/plugins/Vx3DHLTAnalyzer.h"
0011 
0012 #include "FWCore/Utilities/interface/isFinite.h"
0013 #include "FWCore/Framework/interface/LuminosityBlock.h"
0014 
0015 #include <Math/Minimizer.h>
0016 #include <Math/Factory.h>
0017 #include <Math/Functor.h>
0018 
0019 // ### Calling namespaces ###
0020 using namespace std;
0021 using namespace edm;
0022 using namespace reco;
0023 
0024 Vx3DHLTAnalyzer::Vx3DHLTAnalyzer(const ParameterSet& iConfig) {
0025   debugMode = true;
0026   nLumiFit = 2;  // Number of integrated lumis to perform the fit
0027   maxLumiIntegration =
0028       15;  // If failing fits, this is the maximum number of integrated lumis after which a reset is issued
0029   nLumiXaxisRange = 5000;  // Correspond to about 32h of data taking: 32h * 60min * 60s / 23s per lumi-block = 5009
0030   dataFromFit = true;      // The Beam Spot data can be either taken from the histograms or from the fit results
0031   minNentries = 20;        // Minimum number of good vertices to perform the fit
0032   xRange = 0.8;            // [cm]
0033   xStep = 0.001;           // [cm]
0034   yRange = 0.8;            // [cm]
0035   yStep = 0.001;           // [cm]
0036   zRange = 30.;            // [cm]
0037   zStep = 0.04;            // [cm]
0038   VxErrCorr = 1.3;
0039   minVxDoF = 10.;  // Good-vertex selection cut
0040   // For vertex fitter without track-weight: d.o.f. = 2*NTracks - 3
0041   // For vertex fitter with track-weight:    d.o.f. = sum_NTracks(2*track_weight) - 3
0042   minVxWgt = 0.5;  // Good-vertex selection cut
0043   fileName = "BeamPixelResults.txt";
0044 
0045   vertexCollection = consumes<VertexCollection>(
0046       iConfig.getUntrackedParameter<InputTag>("vertexCollection", InputTag("pixelVertices")));
0047   pixelHitCollection = consumes<SiPixelRecHitCollection>(
0048       iConfig.getUntrackedParameter<InputTag>("pixelHitCollection", InputTag("siPixelRecHits")));
0049 
0050   debugMode = iConfig.getParameter<bool>("debugMode");
0051   nLumiFit = iConfig.getParameter<unsigned int>("nLumiFit");
0052   maxLumiIntegration = iConfig.getParameter<unsigned int>("maxLumiIntegration");
0053   nLumiXaxisRange = iConfig.getParameter<unsigned int>("nLumiXaxisRange");
0054   dataFromFit = iConfig.getParameter<bool>("dataFromFit");
0055   minNentries = iConfig.getParameter<unsigned int>("minNentries");
0056   xRange = iConfig.getParameter<double>("xRange");
0057   xStep = iConfig.getParameter<double>("xStep");
0058   yRange = iConfig.getParameter<double>("yRange");
0059   yStep = iConfig.getParameter<double>("yStep");
0060   zRange = iConfig.getParameter<double>("zRange");
0061   zStep = iConfig.getParameter<double>("zStep");
0062   VxErrCorr = iConfig.getParameter<double>("VxErrCorr");
0063   minVxDoF = iConfig.getParameter<double>("minVxDoF");
0064   minVxWgt = iConfig.getParameter<double>("minVxWgt");
0065   fileName = iConfig.getParameter<string>("fileName");
0066 
0067   // ### Set internal variables ###
0068   nParams = 9;  // Number of free parameters in the fit
0069   internalDebug = false;
0070   considerVxCovariance = true;  // Deconvolute vertex covariance matrix
0071   pi = 3.141592653589793238;
0072   // ##############################
0073 }
0074 
0075 Vx3DHLTAnalyzer::~Vx3DHLTAnalyzer() {}
0076 
0077 void Vx3DHLTAnalyzer::analyze(const Event& iEvent, const EventSetup& iSetup) {
0078   Handle<VertexCollection> Vx3DCollection;
0079   iEvent.getByToken(vertexCollection, Vx3DCollection);
0080 
0081   unsigned int i, j;
0082   double det;
0083   VertexType MyVertex;
0084 
0085   if (runNumber != iEvent.id().run()) {
0086     reset("scratch");
0087     runNumber = iEvent.id().run();
0088 
0089     if (debugMode == true) {
0090       stringstream debugFile;
0091       string tmp(fileName);
0092 
0093       if (outputDebugFile.is_open() == true)
0094         outputDebugFile.close();
0095       tmp.erase(strlen(fileName.c_str()) - 4, 4);
0096       debugFile << tmp.c_str() << "_Run" << iEvent.id().run() << ".txt";
0097       outputDebugFile.open(debugFile.str().c_str(), ios::out);
0098       outputDebugFile.close();
0099       outputDebugFile.open(debugFile.str().c_str(), ios::app);
0100     }
0101 
0102     dqmBeginLuminosityBlock(iEvent.getLuminosityBlock(), iSetup);
0103   } else if (beginTimeOfFit != 0) {
0104     totalHits += HitCounter(iEvent);
0105 
0106     if (internalDebug == true) {
0107       cout << "[Vx3DHLTAnalyzer]::\tI found " << totalHits << " pixel hits until now" << endl;
0108       cout << "[Vx3DHLTAnalyzer]::\tIn this event there are " << Vx3DCollection->size() << " vertex cadidates" << endl;
0109     }
0110 
0111     for (vector<Vertex>::const_iterator it3DVx = Vx3DCollection->begin(); it3DVx != Vx3DCollection->end(); it3DVx++) {
0112       if (internalDebug == true) {
0113         cout << "[Vx3DHLTAnalyzer]::\tVertex selections:" << endl;
0114         cout << "[Vx3DHLTAnalyzer]::\tEvent ID = " << iEvent.id() << endl;
0115         cout << "[Vx3DHLTAnalyzer]::\tVertex number = " << it3DVx - Vx3DCollection->begin() << endl;
0116         cout << "[Vx3DHLTAnalyzer]::\tisValid = " << it3DVx->isValid() << endl;
0117         cout << "[Vx3DHLTAnalyzer]::\tisFake = " << it3DVx->isFake() << endl;
0118         cout << "[Vx3DHLTAnalyzer]::\tnodof = " << it3DVx->ndof() << endl;
0119         cout << "[Vx3DHLTAnalyzer]::\ttracksSize = " << it3DVx->tracksSize() << endl;
0120       }
0121 
0122       if ((it3DVx->isValid() == true) && (it3DVx->isFake() == false) && (it3DVx->ndof() >= minVxDoF) &&
0123           (it3DVx->tracksSize() > 0) && ((it3DVx->ndof() + 3.) / ((double)it3DVx->tracksSize()) >= 2. * minVxWgt)) {
0124         for (i = 0; i < DIM; i++) {
0125           for (j = 0; j < DIM; j++) {
0126             MyVertex.Covariance[i][j] = it3DVx->covariance(i, j);
0127             if (isNotFinite(MyVertex.Covariance[i][j]) == true)
0128               break;
0129           }
0130 
0131           if (j != DIM)
0132             break;
0133         }
0134 
0135         if (i == DIM)
0136           det = std::fabs(MyVertex.Covariance[0][0]) *
0137                     (std::fabs(MyVertex.Covariance[1][1]) * std::fabs(MyVertex.Covariance[2][2]) -
0138                      MyVertex.Covariance[1][2] * MyVertex.Covariance[1][2]) -
0139                 MyVertex.Covariance[0][1] * (MyVertex.Covariance[0][1] * std::fabs(MyVertex.Covariance[2][2]) -
0140                                              MyVertex.Covariance[0][2] * MyVertex.Covariance[1][2]) +
0141                 MyVertex.Covariance[0][2] * (MyVertex.Covariance[0][1] * MyVertex.Covariance[1][2] -
0142                                              MyVertex.Covariance[0][2] * std::fabs(MyVertex.Covariance[1][1]));
0143 
0144         if ((i == DIM) && (det > 0.)) {
0145           if (internalDebug == true)
0146             cout << "[Vx3DHLTAnalyzer]::\tVertex accepted !" << endl;
0147 
0148           MyVertex.x = it3DVx->x();
0149           MyVertex.y = it3DVx->y();
0150           MyVertex.z = it3DVx->z();
0151           Vertices.push_back(MyVertex);
0152 
0153           Vx_X->Fill(it3DVx->x());
0154           Vx_Y->Fill(it3DVx->y());
0155           Vx_Z->Fill(it3DVx->z());
0156 
0157           Vx_ZX->Fill(it3DVx->z(), it3DVx->x());
0158           Vx_ZY->Fill(it3DVx->z(), it3DVx->y());
0159           Vx_XY->Fill(it3DVx->x(), it3DVx->y());
0160 
0161           Vx_X_Cum->Fill(it3DVx->x());
0162           Vx_Y_Cum->Fill(it3DVx->y());
0163           Vx_Z_Cum->Fill(it3DVx->z());
0164 
0165           Vx_ZX_Cum->Fill(it3DVx->z(), it3DVx->x());
0166           Vx_ZY_Cum->Fill(it3DVx->z(), it3DVx->y());
0167           Vx_XY_Cum->Fill(it3DVx->x(), it3DVx->y());
0168         } else if (internalDebug == true) {
0169           cout << "[Vx3DHLTAnalyzer]::\tVertex discarded !" << endl;
0170 
0171           for (i = 0; i < DIM; i++)
0172             for (j = 0; j < DIM; j++)
0173               cout << "(i,j) --> " << i << "," << j << " --> " << MyVertex.Covariance[i][j] << endl;
0174         }
0175       } else if (internalDebug == true)
0176         cout << "[Vx3DHLTAnalyzer]::\tVertex discarded !" << endl;
0177     }
0178   }
0179 }
0180 
0181 unsigned int Vx3DHLTAnalyzer::HitCounter(const Event& iEvent) {
0182   Handle<SiPixelRecHitCollection> rechitspixel;
0183   iEvent.getByToken(pixelHitCollection, rechitspixel);
0184 
0185   unsigned int counter = 0;
0186 
0187   for (SiPixelRecHitCollection::const_iterator j = rechitspixel->begin(); j != rechitspixel->end(); j++)
0188     for (edmNew::DetSet<SiPixelRecHit>::const_iterator h = j->begin(); h != j->end(); h++)
0189       counter += h->cluster()->size();
0190 
0191   return counter;
0192 }
0193 
0194 string Vx3DHLTAnalyzer::formatTime(const time_t& t) {
0195   char ts[25];
0196   strftime(ts, sizeof(ts), "%Y.%m.%d %H:%M:%S %Z", gmtime(&t));
0197 
0198   string ts_string(ts);
0199 
0200   return ts_string;
0201 }
0202 
0203 double Vx3DHLTAnalyzer::Gauss3DFunc(const double* par) {
0204   double K[DIM][DIM];  // Covariance Matrix
0205   double M[DIM][DIM];  // K^-1
0206   double det;
0207   double sumlog = 0.;
0208 
0209   //  par[0] = K(0,0) --> Var[X]
0210   //  par[1] = K(1,1) --> Var[Y]
0211   //  par[2] = K(2,2) --> Var[Z]
0212   //  par[3] = K(0,1) = K(1,0) --> Cov[X,Y]
0213   //  par[4] = K(1,2) = K(2,1) --> Cov[Y,Z] --> dy/dz
0214   //  par[5] = K(0,2) = K(2,0) --> Cov[X,Z] --> dx/dz
0215   //  par[6] = mean x
0216   //  par[7] = mean y
0217   //  par[8] = mean z
0218 
0219   counterVx = 0;
0220   for (unsigned int i = 0; i < Vertices.size(); i++) {
0221     if ((std::sqrt((Vertices[i].x - xPos) * (Vertices[i].x - xPos) + (Vertices[i].y - yPos) * (Vertices[i].y - yPos)) <=
0222          maxTransRadius) &&
0223         (std::fabs(Vertices[i].z - zPos) <= maxLongLength)) {
0224       if (considerVxCovariance == true) {
0225         K[0][0] = std::fabs(par[0]) + VxErrCorr * VxErrCorr * std::fabs(Vertices[i].Covariance[0][0]);
0226         K[1][1] = std::fabs(par[1]) + VxErrCorr * VxErrCorr * std::fabs(Vertices[i].Covariance[1][1]);
0227         K[2][2] = std::fabs(par[2]) + VxErrCorr * VxErrCorr * std::fabs(Vertices[i].Covariance[2][2]);
0228         K[0][1] = K[1][0] = par[3] + VxErrCorr * VxErrCorr * Vertices[i].Covariance[0][1];
0229         K[1][2] = K[2][1] = par[4] * (std::fabs(par[2]) - std::fabs(par[1])) - par[5] * par[3] +
0230                             VxErrCorr * VxErrCorr * Vertices[i].Covariance[1][2];
0231         K[0][2] = K[2][0] = par[5] * (std::fabs(par[2]) - std::fabs(par[0])) - par[4] * par[3] +
0232                             VxErrCorr * VxErrCorr * Vertices[i].Covariance[0][2];
0233       } else {
0234         K[0][0] = std::fabs(par[0]);
0235         K[1][1] = std::fabs(par[1]);
0236         K[2][2] = std::fabs(par[2]);
0237         K[0][1] = K[1][0] = par[3];
0238         K[1][2] = K[2][1] = par[4] * (std::fabs(par[2]) - std::fabs(par[1])) - par[5] * par[3];
0239         K[0][2] = K[2][0] = par[5] * (std::fabs(par[2]) - std::fabs(par[0])) - par[4] * par[3];
0240       }
0241 
0242       det = K[0][0] * (K[1][1] * K[2][2] - K[1][2] * K[1][2]) - K[0][1] * (K[0][1] * K[2][2] - K[0][2] * K[1][2]) +
0243             K[0][2] * (K[0][1] * K[1][2] - K[0][2] * K[1][1]);
0244 
0245       M[0][0] = (K[1][1] * K[2][2] - K[1][2] * K[1][2]) / det;
0246       M[1][1] = (K[0][0] * K[2][2] - K[0][2] * K[0][2]) / det;
0247       M[2][2] = (K[0][0] * K[1][1] - K[0][1] * K[0][1]) / det;
0248       M[0][1] = M[1][0] = (K[0][2] * K[1][2] - K[0][1] * K[2][2]) / det;
0249       M[1][2] = M[2][1] = (K[0][2] * K[0][1] - K[1][2] * K[0][0]) / det;
0250       M[0][2] = M[2][0] = (K[0][1] * K[1][2] - K[0][2] * K[1][1]) / det;
0251 
0252       sumlog += double(DIM) * std::log(2. * pi) + std::log(std::fabs(det)) +
0253                 (M[0][0] * (Vertices[i].x - par[6]) * (Vertices[i].x - par[6]) +
0254                  M[1][1] * (Vertices[i].y - par[7]) * (Vertices[i].y - par[7]) +
0255                  M[2][2] * (Vertices[i].z - par[8]) * (Vertices[i].z - par[8]) +
0256                  2. * M[0][1] * (Vertices[i].x - par[6]) * (Vertices[i].y - par[7]) +
0257                  2. * M[1][2] * (Vertices[i].y - par[7]) * (Vertices[i].z - par[8]) +
0258                  2. * M[0][2] * (Vertices[i].x - par[6]) * (Vertices[i].z - par[8]));
0259 
0260       counterVx++;
0261     }
0262   }
0263 
0264   return sumlog;
0265 }
0266 
0267 int Vx3DHLTAnalyzer::MyFit(vector<double>* vals) {
0268   // ############################################
0269   // # RETURN CODE:                             #
0270   // # >0 == NO OK - fit status (MINUIT manual) #
0271   // #  0 == OK                                 #
0272   // # -1 == NO OK - not finite edm             #
0273   // # -2 == NO OK - not enough "minNentries"   #
0274   // # -3 == NO OK - not finite errors          #
0275   // # -4 == NO OK - negative determinant       #
0276   // # -5 == NO OK - maxLumiIntegration reached #
0277   // ############################################
0278 
0279   if ((vals != nullptr) && (vals->size() == nParams * 2)) {
0280     double nSigmaXY = 10.;
0281     double nSigmaZ = 10.;
0282     double parDistanceXY = 1e-3;   // Unit: [cm]
0283     double parDistanceZ = 1e-2;    // Unit: [cm]
0284     double parDistanceddZ = 1e-3;  // Unit: [rad]
0285     double parDistanceCxy = 1e-5;  // Unit: [cm^2]
0286     double bestEdm;
0287 
0288     const unsigned int trials = 4;
0289     double largerDist[trials] = {0.1, 5., 10., 100.};
0290 
0291     double covxz, covyz, det;
0292     double deltaMean;
0293     int bestMovementX = 1;
0294     int bestMovementY = 1;
0295     int bestMovementZ = 1;
0296     int goodData;
0297 
0298     double edm;
0299 
0300     vector<double>::const_iterator it = vals->begin();
0301 
0302     ROOT::Math::Minimizer* Gauss3D = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
0303     Gauss3D->SetErrorDef(1.0);
0304     if (internalDebug == true)
0305       Gauss3D->SetPrintLevel(3);
0306     else
0307       Gauss3D->SetPrintLevel(0);
0308 
0309     ROOT::Math::Functor _Gauss3DFunc(this, &Vx3DHLTAnalyzer::Gauss3DFunc, nParams);
0310     Gauss3D->SetFunction(_Gauss3DFunc);
0311 
0312     if (internalDebug == true)
0313       cout << "[Vx3DHLTAnalyzer]::\t@@@ START FITTING @@@" << endl;
0314 
0315     // @@@ Fit at X-deltaMean | X | X+deltaMean @@@
0316     bestEdm = 1.;
0317     for (int i = 0; i < 3; i++) {
0318       deltaMean = (double(i) - 1.) * std::sqrt(*(it + 0));
0319       if (internalDebug == true)
0320         cout << "[Vx3DHLTAnalyzer]::\tdeltaMean --> " << deltaMean << endl;
0321 
0322       Gauss3D->Clear();
0323 
0324       Gauss3D->SetVariable(0, "var x ", *(it + 0), parDistanceXY * parDistanceXY);
0325       Gauss3D->SetVariable(1, "var y ", *(it + 1), parDistanceXY * parDistanceXY);
0326       Gauss3D->SetVariable(2, "var z ", *(it + 2), parDistanceZ * parDistanceZ);
0327       Gauss3D->SetVariable(3, "cov xy", *(it + 3), parDistanceCxy);
0328       Gauss3D->SetVariable(4, "dydz  ", *(it + 4), parDistanceddZ);
0329       Gauss3D->SetVariable(5, "dxdz  ", *(it + 5), parDistanceddZ);
0330       Gauss3D->SetVariable(6, "mean x", *(it + 6) + deltaMean, parDistanceXY);
0331       Gauss3D->SetVariable(7, "mean y", *(it + 7), parDistanceXY);
0332       Gauss3D->SetVariable(8, "mean z", *(it + 8), parDistanceZ);
0333 
0334       // Set the central positions of the centroid for vertex rejection
0335       xPos = (*vals)[6];
0336       yPos = (*vals)[7];
0337       zPos = (*vals)[8];
0338 
0339       // Set dimensions of the centroid for vertex rejection
0340       maxTransRadius = nSigmaXY * std::sqrt(std::fabs((*vals)[0]) + std::fabs((*vals)[1])) / 2.;
0341       maxLongLength = nSigmaZ * std::sqrt(std::fabs((*vals)[2]));
0342 
0343       Gauss3D->Minimize();
0344       goodData = Gauss3D->Status();
0345       edm = Gauss3D->Edm();
0346 
0347       if (counterVx < minNentries)
0348         goodData = -2;
0349       else if (isNotFinite(edm) == true) {
0350         goodData = -1;
0351         if (internalDebug == true)
0352           cout << "[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl;
0353       } else
0354         for (unsigned int j = 0; j < nParams; j++)
0355           if (isNotFinite(Gauss3D->Errors()[j]) == true) {
0356             goodData = -3;
0357             if (internalDebug == true)
0358               cout << "[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
0359             break;
0360           }
0361       if (goodData == 0) {
0362         covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
0363                 Gauss3D->X()[5] * Gauss3D->X()[3];
0364         covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
0365                 Gauss3D->X()[4] * Gauss3D->X()[3];
0366 
0367         det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
0368               Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
0369               covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
0370         if (det < 0.) {
0371           goodData = -4;
0372           if (internalDebug == true)
0373             cout << "[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl;
0374         }
0375       }
0376 
0377       if ((goodData == 0) && (std::fabs(edm) < bestEdm)) {
0378         bestEdm = edm;
0379         bestMovementX = i;
0380       }
0381     }
0382     if (internalDebug == true)
0383       cout << "[Vx3DHLTAnalyzer]::\tFound bestMovementX --> " << bestMovementX << endl;
0384 
0385     // @@@ Fit at Y-deltaMean | Y | Y+deltaMean @@@
0386     bestEdm = 1.;
0387     for (int i = 0; i < 3; i++) {
0388       deltaMean = (double(i) - 1.) * std::sqrt(*(it + 1));
0389       if (internalDebug == true) {
0390         cout << "[Vx3DHLTAnalyzer]::\tdeltaMean --> " << deltaMean << endl;
0391         cout << "[Vx3DHLTAnalyzer]::\tdeltaMean X --> " << (double(bestMovementX) - 1.) * std::sqrt(*(it + 0)) << endl;
0392       }
0393 
0394       Gauss3D->Clear();
0395 
0396       Gauss3D->SetVariable(0, "var x ", *(it + 0), parDistanceXY * parDistanceXY);
0397       Gauss3D->SetVariable(1, "var y ", *(it + 1), parDistanceXY * parDistanceXY);
0398       Gauss3D->SetVariable(2, "var z ", *(it + 2), parDistanceZ * parDistanceZ);
0399       Gauss3D->SetVariable(3, "cov xy", *(it + 3), parDistanceCxy);
0400       Gauss3D->SetVariable(4, "dydz  ", *(it + 4), parDistanceddZ);
0401       Gauss3D->SetVariable(5, "dxdz  ", *(it + 5), parDistanceddZ);
0402       Gauss3D->SetVariable(6, "mean x", *(it + 6) + (double(bestMovementX) - 1.) * std::sqrt(*(it + 0)), parDistanceXY);
0403       Gauss3D->SetVariable(7, "mean y", *(it + 7) + deltaMean, parDistanceXY);
0404       Gauss3D->SetVariable(8, "mean z", *(it + 8), parDistanceZ);
0405 
0406       // Set the central positions of the centroid for vertex rejection
0407       xPos = Gauss3D->X()[6];
0408       yPos = Gauss3D->X()[7];
0409       zPos = Gauss3D->X()[8];
0410 
0411       // Set dimensions of the centroid for vertex rejection
0412       maxTransRadius = nSigmaXY * std::sqrt(std::fabs(Gauss3D->X()[0]) + std::fabs(Gauss3D->X()[1])) / 2.;
0413       maxLongLength = nSigmaZ * std::sqrt(std::fabs(Gauss3D->X()[2]));
0414 
0415       Gauss3D->Minimize();
0416       goodData = Gauss3D->Status();
0417       edm = Gauss3D->Edm();
0418 
0419       if (counterVx < minNentries)
0420         goodData = -2;
0421       else if (isNotFinite(edm) == true) {
0422         goodData = -1;
0423         if (internalDebug == true)
0424           cout << "[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl;
0425       } else
0426         for (unsigned int j = 0; j < nParams; j++)
0427           if (isNotFinite(Gauss3D->Errors()[j]) == true) {
0428             goodData = -3;
0429             if (internalDebug == true)
0430               cout << "[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
0431             break;
0432           }
0433       if (goodData == 0) {
0434         covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
0435                 Gauss3D->X()[5] * Gauss3D->X()[3];
0436         covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
0437                 Gauss3D->X()[4] * Gauss3D->X()[3];
0438 
0439         det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
0440               Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
0441               covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
0442         if (det < 0.) {
0443           goodData = -4;
0444           if (internalDebug == true)
0445             cout << "[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl;
0446         }
0447       }
0448 
0449       if ((goodData == 0) && (std::fabs(edm) < bestEdm)) {
0450         bestEdm = edm;
0451         bestMovementY = i;
0452       }
0453     }
0454     if (internalDebug == true)
0455       cout << "[Vx3DHLTAnalyzer]::\tFound bestMovementY --> " << bestMovementY << endl;
0456 
0457     // @@@ Fit at Z-deltaMean | Z | Z+deltaMean @@@
0458     bestEdm = 1.;
0459     for (int i = 0; i < 3; i++) {
0460       deltaMean = (double(i) - 1.) * std::sqrt(*(it + 2));
0461       if (internalDebug == true) {
0462         cout << "[Vx3DHLTAnalyzer]::\tdeltaMean --> " << deltaMean << endl;
0463         cout << "[Vx3DHLTAnalyzer]::\tdeltaMean X --> " << (double(bestMovementX) - 1.) * std::sqrt(*(it + 0)) << endl;
0464         cout << "[Vx3DHLTAnalyzer]::\tdeltaMean Y --> " << (double(bestMovementY) - 1.) * std::sqrt(*(it + 1)) << endl;
0465       }
0466 
0467       Gauss3D->Clear();
0468 
0469       Gauss3D->SetVariable(0, "var x ", *(it + 0), parDistanceXY * parDistanceXY);
0470       Gauss3D->SetVariable(1, "var y ", *(it + 1), parDistanceXY * parDistanceXY);
0471       Gauss3D->SetVariable(2, "var z ", *(it + 2), parDistanceZ * parDistanceZ);
0472       Gauss3D->SetVariable(3, "cov xy", *(it + 3), parDistanceCxy);
0473       Gauss3D->SetVariable(4, "dydz  ", *(it + 4), parDistanceddZ);
0474       Gauss3D->SetVariable(5, "dxdz  ", *(it + 5), parDistanceddZ);
0475       Gauss3D->SetVariable(6, "mean x", *(it + 6) + (double(bestMovementX) - 1.) * std::sqrt(*(it + 0)), parDistanceXY);
0476       Gauss3D->SetVariable(7, "mean y", *(it + 7) + (double(bestMovementY) - 1.) * std::sqrt(*(it + 1)), parDistanceXY);
0477       Gauss3D->SetVariable(8, "mean z", *(it + 8) + deltaMean, parDistanceZ);
0478 
0479       // Set the central positions of the centroid for vertex rejection
0480       xPos = Gauss3D->X()[6];
0481       yPos = Gauss3D->X()[7];
0482       zPos = Gauss3D->X()[8];
0483 
0484       // Set dimensions of the centroid for vertex rejection
0485       maxTransRadius = nSigmaXY * std::sqrt(std::fabs(Gauss3D->X()[0]) + std::fabs(Gauss3D->X()[1])) / 2.;
0486       maxLongLength = nSigmaZ * std::sqrt(std::fabs(Gauss3D->X()[2]));
0487 
0488       Gauss3D->Minimize();
0489       goodData = Gauss3D->Status();
0490       edm = Gauss3D->Edm();
0491 
0492       if (counterVx < minNentries)
0493         goodData = -2;
0494       else if (isNotFinite(edm) == true) {
0495         goodData = -1;
0496         if (internalDebug == true)
0497           cout << "[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl;
0498       } else
0499         for (unsigned int j = 0; j < nParams; j++)
0500           if (isNotFinite(Gauss3D->Errors()[j]) == true) {
0501             goodData = -3;
0502             if (internalDebug == true)
0503               cout << "[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
0504             break;
0505           }
0506       if (goodData == 0) {
0507         covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
0508                 Gauss3D->X()[5] * Gauss3D->X()[3];
0509         covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
0510                 Gauss3D->X()[4] * Gauss3D->X()[3];
0511 
0512         det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
0513               Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
0514               covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
0515         if (det < 0.) {
0516           goodData = -4;
0517           if (internalDebug == true)
0518             cout << "[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl;
0519         }
0520       }
0521 
0522       if ((goodData == 0) && (std::fabs(edm) < bestEdm)) {
0523         bestEdm = edm;
0524         bestMovementZ = i;
0525       }
0526     }
0527     if (internalDebug == true)
0528       cout << "[Vx3DHLTAnalyzer]::\tFound bestMovementZ --> " << bestMovementZ << endl;
0529 
0530     Gauss3D->Clear();
0531 
0532     // @@@ FINAL FIT @@@
0533     Gauss3D->SetVariable(0, "var x ", *(it + 0), parDistanceXY * parDistanceXY);
0534     Gauss3D->SetVariable(1, "var y ", *(it + 1), parDistanceXY * parDistanceXY);
0535     Gauss3D->SetVariable(2, "var z ", *(it + 2), parDistanceZ * parDistanceZ);
0536     Gauss3D->SetVariable(3, "cov xy", *(it + 3), parDistanceCxy);
0537     Gauss3D->SetVariable(4, "dydz  ", *(it + 4), parDistanceddZ);
0538     Gauss3D->SetVariable(5, "dxdz  ", *(it + 5), parDistanceddZ);
0539     Gauss3D->SetVariable(6, "mean x", *(it + 6) + (double(bestMovementX) - 1.) * std::sqrt(*(it + 0)), parDistanceXY);
0540     Gauss3D->SetVariable(7, "mean y", *(it + 7) + (double(bestMovementY) - 1.) * std::sqrt(*(it + 1)), parDistanceXY);
0541     Gauss3D->SetVariable(8, "mean z", *(it + 8) + (double(bestMovementZ) - 1.) * std::sqrt(*(it + 2)), parDistanceZ);
0542 
0543     // Set the central positions of the centroid for vertex rejection
0544     xPos = (*vals)[6];
0545     yPos = (*vals)[7];
0546     zPos = (*vals)[8];
0547 
0548     // Set dimensions of the centroid for vertex rejection
0549     maxTransRadius = nSigmaXY * std::sqrt(std::fabs((*vals)[0]) + std::fabs((*vals)[1])) / 2.;
0550     maxLongLength = nSigmaZ * std::sqrt(std::fabs((*vals)[2]));
0551 
0552     Gauss3D->Minimize();
0553     goodData = Gauss3D->Status();
0554     edm = Gauss3D->Edm();
0555 
0556     if (counterVx < minNentries)
0557       goodData = -2;
0558     else if (isNotFinite(edm) == true) {
0559       goodData = -1;
0560       if (internalDebug == true)
0561         cout << "[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl;
0562     } else
0563       for (unsigned int j = 0; j < nParams; j++)
0564         if (isNotFinite(Gauss3D->Errors()[j]) == true) {
0565           goodData = -3;
0566           if (internalDebug == true)
0567             cout << "[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
0568           break;
0569         }
0570     if (goodData == 0) {
0571       covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
0572               Gauss3D->X()[5] * Gauss3D->X()[3];
0573       covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
0574               Gauss3D->X()[4] * Gauss3D->X()[3];
0575 
0576       det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
0577             Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
0578             covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
0579       if (det < 0.) {
0580         goodData = -4;
0581         if (internalDebug == true)
0582           cout << "[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl;
0583       }
0584     }
0585 
0586     // @@@ FIT WITH DIFFERENT PARAMETER DISTANCES @@@
0587     for (unsigned int i = 0; i < trials; i++) {
0588       if ((goodData != 0) && (goodData != -2)) {
0589         Gauss3D->Clear();
0590 
0591         if (internalDebug == true)
0592           cout << "[Vx3DHLTAnalyzer]::\tFIT WITH DIFFERENT PARAMETER DISTANCES - STEP " << i + 1 << endl;
0593 
0594         Gauss3D->SetVariable(0, "var x ", *(it + 0), parDistanceXY * parDistanceXY * largerDist[i]);
0595         Gauss3D->SetVariable(1, "var y ", *(it + 1), parDistanceXY * parDistanceXY * largerDist[i]);
0596         Gauss3D->SetVariable(2, "var z ", *(it + 2), parDistanceZ * parDistanceZ * largerDist[i]);
0597         Gauss3D->SetVariable(3, "cov xy", *(it + 3), parDistanceCxy * largerDist[i]);
0598         Gauss3D->SetVariable(4, "dydz  ", *(it + 4), parDistanceddZ * largerDist[i]);
0599         Gauss3D->SetVariable(5, "dxdz  ", *(it + 5), parDistanceddZ * largerDist[i]);
0600         Gauss3D->SetVariable(6,
0601                              "mean x",
0602                              *(it + 6) + (double(bestMovementX) - 1.) * std::sqrt(*(it + 0)),
0603                              parDistanceXY * largerDist[i]);
0604         Gauss3D->SetVariable(7,
0605                              "mean y",
0606                              *(it + 7) + (double(bestMovementY) - 1.) * std::sqrt(*(it + 1)),
0607                              parDistanceXY * largerDist[i]);
0608         Gauss3D->SetVariable(
0609             8, "mean z", *(it + 8) + (double(bestMovementZ) - 1.) * std::sqrt(*(it + 2)), parDistanceZ * largerDist[i]);
0610 
0611         // Set the central positions of the centroid for vertex rejection
0612         xPos = Gauss3D->X()[6];
0613         yPos = Gauss3D->X()[7];
0614         zPos = Gauss3D->X()[8];
0615 
0616         // Set dimensions of the centroid for vertex rejection
0617         maxTransRadius = nSigmaXY * std::sqrt(std::fabs(Gauss3D->X()[0]) + std::fabs(Gauss3D->X()[1])) / 2.;
0618         maxLongLength = nSigmaZ * std::sqrt(std::fabs(Gauss3D->X()[2]));
0619 
0620         Gauss3D->Minimize();
0621         goodData = Gauss3D->Status();
0622         edm = Gauss3D->Edm();
0623 
0624         if (counterVx < minNentries)
0625           goodData = -2;
0626         else if (isNotFinite(edm) == true) {
0627           goodData = -1;
0628           if (internalDebug == true)
0629             cout << "[Vx3DHLTAnalyzer]::\tNot finite edm !" << endl;
0630         } else
0631           for (unsigned int j = 0; j < nParams; j++)
0632             if (isNotFinite(Gauss3D->Errors()[j]) == true) {
0633               goodData = -3;
0634               if (internalDebug == true)
0635                 cout << "[Vx3DHLTAnalyzer]::\tNot finite errors !" << endl;
0636               break;
0637             }
0638         if (goodData == 0) {
0639           covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
0640                   Gauss3D->X()[5] * Gauss3D->X()[3];
0641           covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
0642                   Gauss3D->X()[4] * Gauss3D->X()[3];
0643 
0644           det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
0645                 Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
0646                 covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
0647           if (det < 0.) {
0648             goodData = -4;
0649             if (internalDebug == true)
0650               cout << "[Vx3DHLTAnalyzer]::\tNegative determinant !" << endl;
0651           }
0652         }
0653       } else
0654         break;
0655     }
0656 
0657     if (goodData == 0)
0658       for (unsigned int i = 0; i < nParams; i++) {
0659         vals->operator[](i) = Gauss3D->X()[i];
0660         vals->operator[](i + nParams) = Gauss3D->Errors()[i];
0661       }
0662 
0663     delete Gauss3D;
0664     return goodData;
0665   }
0666 
0667   return -1;
0668 }
0669 
0670 void Vx3DHLTAnalyzer::reset(string ResetType) {
0671   if ((debugMode == true) && (outputDebugFile.is_open() == true)) {
0672     outputDebugFile << "Runnumber " << runNumber << endl;
0673     outputDebugFile << "BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) << " " << (beginTimeOfFit >> 32) << endl;
0674     outputDebugFile << "BeginLumiRange " << beginLumiOfFit << endl;
0675     outputDebugFile << "EndTimeOfFit " << formatTime(endTimeOfFit >> 32) << " " << (endTimeOfFit >> 32) << endl;
0676     outputDebugFile << "EndLumiRange " << endLumiOfFit << endl;
0677     outputDebugFile << "LumiCounter " << lumiCounter << endl;
0678     outputDebugFile << "LastLumiOfFit " << lastLumiOfFit << endl;
0679   }
0680 
0681   if (ResetType == "scratch") {
0682     runNumber = 0;
0683     numberGoodFits = 0;
0684     numberFits = 0;
0685     lastLumiOfFit = 0;
0686 
0687     Vx_X->Reset();
0688     Vx_Y->Reset();
0689     Vx_Z->Reset();
0690 
0691     Vx_X_Fit->Reset();
0692     Vx_Y_Fit->Reset();
0693     Vx_Z_Fit->Reset();
0694 
0695     Vx_ZX->Reset();
0696     Vx_ZY->Reset();
0697     Vx_XY->Reset();
0698 
0699     Vx_X_Cum->Reset();
0700     Vx_Y_Cum->Reset();
0701     Vx_Z_Cum->Reset();
0702 
0703     Vx_ZX_Cum->Reset();
0704     Vx_ZY_Cum->Reset();
0705     Vx_XY_Cum->Reset();
0706 
0707     mXlumi->Reset();
0708     mYlumi->Reset();
0709     mZlumi->Reset();
0710 
0711     sXlumi->Reset();
0712     sYlumi->Reset();
0713     sZlumi->Reset();
0714 
0715     dxdzlumi->Reset();
0716     dydzlumi->Reset();
0717 
0718     hitCounter->Reset();
0719     goodVxCounter->Reset();
0720     statusCounter->Reset();
0721     fitResults->Reset();
0722 
0723     reportSummary->Fill(-1);
0724     reportSummaryMap->getTH1()->SetBinContent(1, 1, -1);
0725 
0726     Vertices.clear();
0727 
0728     lumiCounter = 0;
0729     totalHits = 0;
0730     beginTimeOfFit = 0;
0731     endTimeOfFit = 0;
0732     beginLumiOfFit = 0;
0733     endLumiOfFit = 0;
0734 
0735     if (internalDebug == true)
0736       cout << "[Vx3DHLTAnalyzer]::\tReset issued: scratch" << endl;
0737     if ((debugMode == true) && (outputDebugFile.is_open() == true))
0738       outputDebugFile << "Reset -scratch- issued\n" << endl;
0739   } else if (ResetType == "whole") {
0740     Vx_X->Reset();
0741     Vx_Y->Reset();
0742     Vx_Z->Reset();
0743 
0744     Vx_ZX->Reset();
0745     Vx_ZY->Reset();
0746     Vx_XY->Reset();
0747 
0748     Vertices.clear();
0749 
0750     lumiCounter = 0;
0751     totalHits = 0;
0752     beginTimeOfFit = 0;
0753     endTimeOfFit = 0;
0754     beginLumiOfFit = 0;
0755     endLumiOfFit = 0;
0756 
0757     if (internalDebug == true)
0758       cout << "[Vx3DHLTAnalyzer]::\tReset issued: whole" << endl;
0759     if ((debugMode == true) && (outputDebugFile.is_open() == true))
0760       outputDebugFile << "Reset -whole- issued\n" << endl;
0761   } else if (ResetType == "fit") {
0762     Vx_X_Fit->Reset();
0763     Vx_Y_Fit->Reset();
0764     Vx_Z_Fit->Reset();
0765 
0766     if (internalDebug == true)
0767       cout << "[Vx3DHLTAnalyzer]::\tReset issued: fit" << endl;
0768     if ((debugMode == true) && (outputDebugFile.is_open() == true))
0769       outputDebugFile << "Reset -fit- issued\n" << endl;
0770   } else if (ResetType == "hitCounter") {
0771     totalHits = 0;
0772 
0773     if (internalDebug == true)
0774       cout << "[Vx3DHLTAnalyzer]::\tReset issued: hitCounter" << endl;
0775     if ((debugMode == true) && (outputDebugFile.is_open() == true))
0776       outputDebugFile << "Reset -hitCounter- issued\n" << endl;
0777   }
0778 }
0779 
0780 void Vx3DHLTAnalyzer::writeToFile(vector<double>* vals,
0781                                   TimeValue_t BeginTimeOfFit,
0782                                   TimeValue_t EndTimeOfFit,
0783                                   unsigned int BeginLumiOfFit,
0784                                   unsigned int EndLumiOfFit,
0785                                   int dataType) {
0786   stringstream BufferString;
0787   BufferString.precision(5);
0788 
0789   outputFile.open(fileName.c_str(), ios::out);
0790 
0791   if ((outputFile.is_open() == true) && (vals != nullptr) && (vals->size() == (nParams - 1) * 2)) {
0792     vector<double>::const_iterator it = vals->begin();
0793 
0794     outputFile << "Runnumber " << runNumber << endl;
0795     outputFile << "BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) << " " << (beginTimeOfFit >> 32) << endl;
0796     outputFile << "EndTimeOfFit " << formatTime(endTimeOfFit >> 32) << " " << (endTimeOfFit >> 32) << endl;
0797     outputFile << "LumiRange " << beginLumiOfFit << " - " << endLumiOfFit << endl;
0798     outputFile << "Type " << dataType << endl;
0799     // 3D Vertexing with Pixel Tracks:
0800     // Good data = Type  3
0801     // Bad data  = Type -1
0802 
0803     BufferString << *(it + 0);
0804     outputFile << "X0 " << BufferString.str().c_str() << endl;
0805     BufferString.str("");
0806 
0807     BufferString << *(it + 1);
0808     outputFile << "Y0 " << BufferString.str().c_str() << endl;
0809     BufferString.str("");
0810 
0811     BufferString << *(it + 2);
0812     outputFile << "Z0 " << BufferString.str().c_str() << endl;
0813     BufferString.str("");
0814 
0815     BufferString << *(it + 5);
0816     outputFile << "sigmaZ0 " << BufferString.str().c_str() << endl;
0817     BufferString.str("");
0818 
0819     BufferString << *(it + 6);
0820     outputFile << "dxdz " << BufferString.str().c_str() << endl;
0821     BufferString.str("");
0822 
0823     BufferString << *(it + 7);
0824     outputFile << "dydz " << BufferString.str().c_str() << endl;
0825     BufferString.str("");
0826 
0827     BufferString << *(it + 3);
0828     outputFile << "BeamWidthX " << BufferString.str().c_str() << endl;
0829     BufferString.str("");
0830 
0831     BufferString << *(it + 4);
0832     outputFile << "BeamWidthY " << BufferString.str().c_str() << endl;
0833     BufferString.str("");
0834 
0835     outputFile << "Cov(0,j) " << *(it + 8) << " 0 0 0 0 0 0" << endl;
0836     outputFile << "Cov(1,j) 0 " << *(it + 9) << " 0 0 0 0 0" << endl;
0837     outputFile << "Cov(2,j) 0 0 " << *(it + 10) << " 0 0 0 0" << endl;
0838     outputFile << "Cov(3,j) 0 0 0 " << *(it + 13) << " 0 0 0" << endl;
0839     outputFile << "Cov(4,j) 0 0 0 0 " << *(it + 14) << " 0 0" << endl;
0840     outputFile << "Cov(5,j) 0 0 0 0 0 " << *(it + 15) << " 0" << endl;
0841     outputFile << "Cov(6,j) 0 0 0 0 0 0 "
0842                << ((*(it + 11)) + (*(it + 12)) + 2. * std::sqrt((*(it + 11)) * (*(it + 12)))) / 4. << endl;
0843 
0844     outputFile << "EmittanceX 0" << endl;
0845     outputFile << "EmittanceY 0" << endl;
0846     outputFile << "BetaStar 0" << endl;
0847   }
0848   outputFile.close();
0849 
0850   if ((debugMode == true) && (outputDebugFile.is_open() == true) && (vals != nullptr) &&
0851       (vals->size() == (nParams - 1) * 2)) {
0852     vector<double>::const_iterator it = vals->begin();
0853 
0854     outputDebugFile << "Runnumber " << runNumber << endl;
0855     outputDebugFile << "BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) << " " << (beginTimeOfFit >> 32) << endl;
0856     outputDebugFile << "EndTimeOfFit " << formatTime(endTimeOfFit >> 32) << " " << (endTimeOfFit >> 32) << endl;
0857     outputDebugFile << "LumiRange " << beginLumiOfFit << " - " << endLumiOfFit << endl;
0858     outputDebugFile << "Type " << dataType << endl;
0859     // 3D Vertexing with Pixel Tracks:
0860     // Good data = Type  3
0861     // Bad data  = Type -1
0862 
0863     BufferString << *(it + 0);
0864     outputDebugFile << "X0 " << BufferString.str().c_str() << endl;
0865     BufferString.str("");
0866 
0867     BufferString << *(it + 1);
0868     outputDebugFile << "Y0 " << BufferString.str().c_str() << endl;
0869     BufferString.str("");
0870 
0871     BufferString << *(it + 2);
0872     outputDebugFile << "Z0 " << BufferString.str().c_str() << endl;
0873     BufferString.str("");
0874 
0875     BufferString << *(it + 5);
0876     outputDebugFile << "sigmaZ0 " << BufferString.str().c_str() << endl;
0877     BufferString.str("");
0878 
0879     BufferString << *(it + 6);
0880     outputDebugFile << "dxdz " << BufferString.str().c_str() << endl;
0881     BufferString.str("");
0882 
0883     BufferString << *(it + 7);
0884     outputDebugFile << "dydz " << BufferString.str().c_str() << endl;
0885     BufferString.str("");
0886 
0887     BufferString << *(it + 3);
0888     outputDebugFile << "BeamWidthX " << BufferString.str().c_str() << endl;
0889     BufferString.str("");
0890 
0891     BufferString << *(it + 4);
0892     outputDebugFile << "BeamWidthY " << BufferString.str().c_str() << endl;
0893     BufferString.str("");
0894 
0895     outputDebugFile << "Cov(0,j) " << *(it + 8) << " 0 0 0 0 0 0" << endl;
0896     outputDebugFile << "Cov(1,j) 0 " << *(it + 9) << " 0 0 0 0 0" << endl;
0897     outputDebugFile << "Cov(2,j) 0 0 " << *(it + 10) << " 0 0 0 0" << endl;
0898     outputDebugFile << "Cov(3,j) 0 0 0 " << *(it + 13) << " 0 0 0" << endl;
0899     outputDebugFile << "Cov(4,j) 0 0 0 0 " << *(it + 14) << " 0 0" << endl;
0900     outputDebugFile << "Cov(5,j) 0 0 0 0 0 " << *(it + 15) << " 0" << endl;
0901     outputDebugFile << "Cov(6,j) 0 0 0 0 0 0 "
0902                     << ((*(it + 11)) + (*(it + 12)) + 2. * std::sqrt((*(it + 11)) * (*(it + 12)))) / 4. << endl;
0903 
0904     outputDebugFile << "EmittanceX 0" << endl;
0905     outputDebugFile << "EmittanceY 0" << endl;
0906     outputDebugFile << "BetaStar 0" << endl;
0907 
0908     outputDebugFile << "\nUsed vertices: " << counterVx << "\n" << endl;
0909   }
0910 }
0911 
0912 void Vx3DHLTAnalyzer::printFitParams(const vector<double>& fitResults) {
0913   cout << "var x -->  " << fitResults[0] << " +/- " << fitResults[0 + nParams] << endl;
0914   cout << "var y -->  " << fitResults[1] << " +/- " << fitResults[1 + nParams] << endl;
0915   cout << "var z -->  " << fitResults[2] << " +/- " << fitResults[2 + nParams] << endl;
0916   cout << "cov xy --> " << fitResults[3] << " +/- " << fitResults[3 + nParams] << endl;
0917   cout << "dydz   --> " << fitResults[4] << " +/- " << fitResults[4 + nParams] << endl;
0918   cout << "dxdz   --> " << fitResults[5] << " +/- " << fitResults[5 + nParams] << endl;
0919   cout << "mean x --> " << fitResults[6] << " +/- " << fitResults[6 + nParams] << endl;
0920   cout << "mean y --> " << fitResults[7] << " +/- " << fitResults[7 + nParams] << endl;
0921   cout << "mean z --> " << fitResults[8] << " +/- " << fitResults[8 + nParams] << endl;
0922 }
0923 
0924 void Vx3DHLTAnalyzer::dqmBeginLuminosityBlock(const LuminosityBlock& lumiBlock, const EventSetup& iSetup) {
0925   // @@@ If statement to avoid problems with non-sequential lumisections @@@
0926   if ((lumiCounter == 0) && (lumiBlock.luminosityBlock() > lastLumiOfFit)) {
0927     beginTimeOfFit = lumiBlock.beginTime().value();
0928     beginLumiOfFit = lumiBlock.luminosityBlock();
0929     lumiCounter++;
0930   } else if ((lumiCounter != 0) && (lumiBlock.luminosityBlock() >= (beginLumiOfFit + lumiCounter)))
0931     lumiCounter++;
0932   else
0933     reset("scratch");
0934 }
0935 
0936 void Vx3DHLTAnalyzer::dqmEndLuminosityBlock(const LuminosityBlock& lumiBlock, const EventSetup& iSetup) {
0937   stringstream histTitle;
0938   double minXfit, maxXfit;
0939   int goodData;
0940 
0941   if ((nLumiFit != 0) && (lumiCounter % nLumiFit == 0) && (beginTimeOfFit != 0) && (runNumber != 0)) {
0942     endTimeOfFit = lumiBlock.endTime().value();
0943     endLumiOfFit = lumiBlock.luminosityBlock();
0944     lastLumiOfFit = endLumiOfFit;
0945     vector<double> vals;
0946 
0947     hitCounter->getTH1()->SetBinContent(lastLumiOfFit, (double)totalHits);
0948     hitCounter->getTH1()->SetBinError(
0949         lastLumiOfFit,
0950         (totalHits != 0 ? 1.
0951                         : 0.));  // It's not sqrt(n) because we want to weight all entries in the same way for the fit
0952 
0953     if (dataFromFit == true) {
0954       vector<double> fitResults;
0955 
0956       fitResults.push_back(Vx_X->getTH1()->GetRMS() * Vx_X->getTH1()->GetRMS());
0957       fitResults.push_back(Vx_Y->getTH1()->GetRMS() * Vx_Y->getTH1()->GetRMS());
0958       fitResults.push_back(Vx_Z->getTH1()->GetRMS() * Vx_Z->getTH1()->GetRMS());
0959       fitResults.push_back(0.0);
0960       fitResults.push_back(0.0);
0961       fitResults.push_back(0.0);
0962       fitResults.push_back(Vx_X->getTH1()->GetMean());
0963       fitResults.push_back(Vx_Y->getTH1()->GetMean());
0964       fitResults.push_back(Vx_Z->getTH1()->GetMean());
0965       for (unsigned int i = 0; i < nParams; i++)
0966         fitResults.push_back(0.0);
0967 
0968       if (internalDebug == true) {
0969         cout << "[Vx3DHLTAnalyzer]::\t@@@ Beam Spot parameters - prefit @@@" << endl;
0970 
0971         printFitParams(fitResults);
0972 
0973         cout << "Runnumber " << runNumber << endl;
0974         cout << "BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) << " " << (beginTimeOfFit >> 32) << endl;
0975         cout << "EndTimeOfFit " << formatTime(endTimeOfFit >> 32) << " " << (endTimeOfFit >> 32) << endl;
0976         cout << "LumiRange " << beginLumiOfFit << " - " << endLumiOfFit << endl;
0977       }
0978 
0979       goodData = MyFit(&fitResults);
0980 
0981       if (internalDebug == true) {
0982         cout << "[Vx3DHLTAnalyzer]::\t@@@ Beam Spot parameters - postfit @@@" << endl;
0983 
0984         printFitParams(fitResults);
0985 
0986         cout << "goodData --> " << goodData << endl;
0987         cout << "Used vertices --> " << counterVx << endl;
0988       }
0989 
0990       if (goodData == 0) {
0991         vals.push_back(fitResults[6]);
0992         vals.push_back(fitResults[7]);
0993         vals.push_back(fitResults[8]);
0994         vals.push_back(std::sqrt(std::fabs(fitResults[0])));
0995         vals.push_back(std::sqrt(std::fabs(fitResults[1])));
0996         vals.push_back(std::sqrt(std::fabs(fitResults[2])));
0997         vals.push_back(fitResults[5]);
0998         vals.push_back(fitResults[4]);
0999 
1000         vals.push_back(std::pow(fitResults[6 + nParams], 2.));
1001         vals.push_back(std::pow(fitResults[7 + nParams], 2.));
1002         vals.push_back(std::pow(fitResults[8 + nParams], 2.));
1003         vals.push_back(std::pow(std::fabs(fitResults[0 + nParams]) / (2. * std::sqrt(std::fabs(fitResults[0]))), 2.));
1004         vals.push_back(std::pow(std::fabs(fitResults[1 + nParams]) / (2. * std::sqrt(std::fabs(fitResults[1]))), 2.));
1005         vals.push_back(std::pow(std::fabs(fitResults[2 + nParams]) / (2. * std::sqrt(std::fabs(fitResults[2]))), 2.));
1006         vals.push_back(std::pow(fitResults[5 + nParams], 2.));
1007         vals.push_back(std::pow(fitResults[4 + nParams], 2.));
1008       } else
1009         for (unsigned int i = 0; i < (nParams - 1) * 2; i++)
1010           vals.push_back(0.0);
1011 
1012       fitResults.clear();
1013     } else {
1014       counterVx = Vx_X->getTH1F()->GetEntries();
1015 
1016       if (Vx_X->getTH1F()->GetEntries() >= minNentries) {
1017         goodData = 0;
1018 
1019         vals.push_back(Vx_X->getTH1F()->GetMean());
1020         vals.push_back(Vx_Y->getTH1F()->GetMean());
1021         vals.push_back(Vx_Z->getTH1F()->GetMean());
1022         vals.push_back(Vx_X->getTH1F()->GetRMS());
1023         vals.push_back(Vx_Y->getTH1F()->GetRMS());
1024         vals.push_back(Vx_Z->getTH1F()->GetRMS());
1025         vals.push_back(0.0);
1026         vals.push_back(0.0);
1027 
1028         vals.push_back(std::pow(Vx_X->getTH1F()->GetMeanError(), 2.));
1029         vals.push_back(std::pow(Vx_Y->getTH1F()->GetMeanError(), 2.));
1030         vals.push_back(std::pow(Vx_Z->getTH1F()->GetMeanError(), 2.));
1031         vals.push_back(std::pow(Vx_X->getTH1F()->GetRMSError(), 2.));
1032         vals.push_back(std::pow(Vx_Y->getTH1F()->GetRMSError(), 2.));
1033         vals.push_back(std::pow(Vx_Z->getTH1F()->GetRMSError(), 2.));
1034         vals.push_back(0.0);
1035         vals.push_back(0.0);
1036       } else {
1037         goodData = -2;
1038         for (unsigned int i = 0; i < (nParams - 1) * 2; i++)
1039           vals.push_back(0.0);
1040       }
1041     }
1042 
1043     // vals[0]  = X0
1044     // vals[1]  = Y0
1045     // vals[2]  = Z0
1046     // vals[3]  = sigmaX0
1047     // vals[4]  = sigmaY0
1048     // vals[5]  = sigmaZ0
1049     // vals[6]  = dxdz
1050     // vals[7]  = dydz
1051 
1052     // vals[8]  = err^2 X0
1053     // vals[9]  = err^2 Y0
1054     // vals[10] = err^2 Z0
1055     // vals[11] = err^2 sigmaX0
1056     // vals[12] = err^2 sigmaY0
1057     // vals[13] = err^2 sigmaZ0
1058     // vals[14] = err^2 dxdz
1059     // vals[15] = err^2 dydz
1060 
1061     numberFits++;
1062     writeToFile(&vals, beginTimeOfFit, endTimeOfFit, beginLumiOfFit, endLumiOfFit, 3);
1063     if (internalDebug == true)
1064       cout << "[Vx3DHLTAnalyzer]::\tUsed vertices: " << counterVx << endl;
1065 
1066     statusCounter->getTH1()->SetBinContent(lastLumiOfFit, (double)goodData);
1067     statusCounter->getTH1()->SetBinError(lastLumiOfFit, 1e-3);
1068 
1069     // Copy vertex position histograms into to-fit histograms
1070     if (goodData == 0)
1071       reset("fit");
1072     else if (lumiCounter >= maxLumiIntegration) {
1073       reset("fit");
1074       reset("whole");
1075     }
1076 
1077     for (int i = 0; i < Vx_X_Fit->getTH1()->GetNbinsX(); i++) {
1078       Vx_X_Fit->getTH1()->SetBinContent(
1079           i + 1, Vx_X_Fit->getTH1()->GetBinContent(i + 1) + Vx_X->getTH1()->GetBinContent(i + 1));
1080       Vx_X_Fit->getTH1()->SetBinError(i + 1, sqrt(Vx_X_Fit->getTH1()->GetBinContent(i + 1)));
1081     }
1082 
1083     for (int i = 0; i < Vx_Y_Fit->getTH1()->GetNbinsX(); i++) {
1084       Vx_Y_Fit->getTH1()->SetBinContent(
1085           i + 1, Vx_Y_Fit->getTH1()->GetBinContent(i + 1) + Vx_Y->getTH1()->GetBinContent(i + 1));
1086       Vx_Y_Fit->getTH1()->SetBinError(i + 1, sqrt(Vx_Y_Fit->getTH1()->GetBinContent(i + 1)));
1087     }
1088 
1089     for (int i = 0; i < Vx_Z_Fit->getTH1()->GetNbinsX(); i++) {
1090       Vx_Z_Fit->getTH1()->SetBinContent(
1091           i + 1, Vx_Z_Fit->getTH1()->GetBinContent(i + 1) + Vx_Z->getTH1()->GetBinContent(i + 1));
1092       Vx_Z_Fit->getTH1()->SetBinError(i + 1, sqrt(Vx_Z_Fit->getTH1()->GetBinContent(i + 1)));
1093     }
1094 
1095     // Check data quality
1096     if (goodData == 0) {
1097       numberGoodFits++;
1098 
1099       histTitle << "Ongoing: fitted lumis " << beginLumiOfFit << " - " << endLumiOfFit;
1100       reset("whole");
1101     } else {
1102       if (goodData == -2)
1103         histTitle << "Ongoing: not enough evts (" << lumiCounter << " - " << maxLumiIntegration << " lumis)";
1104       else
1105         histTitle << "Ongoing: temporary problems (" << lumiCounter << " - " << maxLumiIntegration << " lumis)";
1106 
1107       if (lumiCounter >= maxLumiIntegration) {
1108         statusCounter->getTH1()->SetBinContent(lastLumiOfFit, -5);
1109         statusCounter->getTH1()->SetBinError(lastLumiOfFit, 1e-3);
1110       } else
1111         reset("hitCounter");
1112     }
1113 
1114     reportSummary->Fill((numberFits != 0 ? ((double)numberGoodFits) / ((double)numberFits) : -1));
1115     reportSummaryMap->getTH1()->SetBinContent(
1116         1, 1, (numberFits != 0 ? ((double)numberGoodFits) / ((double)numberFits) : -1));
1117 
1118     fitResults->setAxisTitle(histTitle.str(), 1);
1119 
1120     fitResults->setBinContent(1, 9, vals[0]);
1121     fitResults->setBinContent(1, 8, vals[1]);
1122     fitResults->setBinContent(1, 7, vals[2]);
1123     fitResults->setBinContent(1, 6, vals[3]);
1124     fitResults->setBinContent(1, 5, vals[4]);
1125     fitResults->setBinContent(1, 4, vals[5]);
1126     fitResults->setBinContent(1, 3, vals[6]);
1127     fitResults->setBinContent(1, 2, vals[7]);
1128     fitResults->setBinContent(1, 1, counterVx);
1129 
1130     fitResults->setBinContent(2, 9, std::sqrt(vals[8]));
1131     fitResults->setBinContent(2, 8, std::sqrt(vals[9]));
1132     fitResults->setBinContent(2, 7, std::sqrt(vals[10]));
1133     fitResults->setBinContent(2, 6, std::sqrt(vals[11]));
1134     fitResults->setBinContent(2, 5, std::sqrt(vals[12]));
1135     fitResults->setBinContent(2, 4, std::sqrt(vals[13]));
1136     fitResults->setBinContent(2, 3, std::sqrt(vals[14]));
1137     fitResults->setBinContent(2, 2, std::sqrt(vals[15]));
1138     fitResults->setBinContent(2, 1, std::sqrt(counterVx));
1139 
1140     // Linear fit to the historical plots
1141     TF1* myLinFit = new TF1(
1142         "myLinFit", "[0] + [1]*x", mXlumi->getTH1()->GetXaxis()->GetXmin(), mXlumi->getTH1()->GetXaxis()->GetXmax());
1143     myLinFit->SetLineColor(2);
1144     myLinFit->SetLineWidth(2);
1145     myLinFit->SetParName(0, "Inter.");
1146     myLinFit->SetParName(1, "Slope");
1147 
1148     mXlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[0]);
1149     mXlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[8]));
1150     myLinFit->SetParameter(0, mXlumi->getTH1()->GetMean(2));
1151     myLinFit->SetParameter(1, 0.0);
1152     mXlumi->getTH1()->Fit(myLinFit, "QR");
1153 
1154     mYlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[1]);
1155     mYlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[9]));
1156     myLinFit->SetParameter(0, mYlumi->getTH1()->GetMean(2));
1157     myLinFit->SetParameter(1, 0.0);
1158     mYlumi->getTH1()->Fit(myLinFit, "QR");
1159 
1160     mZlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[2]);
1161     mZlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[10]));
1162     myLinFit->SetParameter(0, mZlumi->getTH1()->GetMean(2));
1163     myLinFit->SetParameter(1, 0.0);
1164     mZlumi->getTH1()->Fit(myLinFit, "QR");
1165 
1166     sXlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[3]);
1167     sXlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[11]));
1168     myLinFit->SetParameter(0, sXlumi->getTH1()->GetMean(2));
1169     myLinFit->SetParameter(1, 0.0);
1170     sXlumi->getTH1()->Fit(myLinFit, "QR");
1171 
1172     sYlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[4]);
1173     sYlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[12]));
1174     myLinFit->SetParameter(0, sYlumi->getTH1()->GetMean(2));
1175     myLinFit->SetParameter(1, 0.0);
1176     sYlumi->getTH1()->Fit(myLinFit, "QR");
1177 
1178     sZlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[5]);
1179     sZlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[13]));
1180     myLinFit->SetParameter(0, sZlumi->getTH1()->GetMean(2));
1181     myLinFit->SetParameter(1, 0.0);
1182     sZlumi->getTH1()->Fit(myLinFit, "QR");
1183 
1184     dxdzlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[6]);
1185     dxdzlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[14]));
1186     myLinFit->SetParameter(0, dxdzlumi->getTH1()->GetMean(2));
1187     myLinFit->SetParameter(1, 0.0);
1188     dxdzlumi->getTH1()->Fit(myLinFit, "QR");
1189 
1190     dydzlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[7]);
1191     dydzlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[15]));
1192     myLinFit->SetParameter(0, dydzlumi->getTH1()->GetMean(2));
1193     myLinFit->SetParameter(1, 0.0);
1194     dydzlumi->getTH1()->Fit(myLinFit, "QR");
1195 
1196     myLinFit->SetParameter(0, hitCounter->getTH1()->GetMean(2));
1197     myLinFit->SetParameter(1, 0.0);
1198     hitCounter->getTH1()->Fit(myLinFit, "QR");
1199 
1200     goodVxCounter->getTH1()->SetBinContent(lastLumiOfFit, (double)counterVx);
1201     goodVxCounter->getTH1()->SetBinError(
1202         lastLumiOfFit,
1203         (counterVx != 0 ? 1.
1204                         : 0.));  // It's not sqrt(n) because we want to weight all entries in the same way for the fit
1205     myLinFit->SetParameter(0, goodVxCounter->getTH1()->GetMean(2));
1206     myLinFit->SetParameter(1, 0.0);
1207     goodVxCounter->getTH1()->Fit(myLinFit, "QR");
1208 
1209     delete myLinFit;
1210     vals.clear();
1211 
1212     // Gaussian fit to 1D vertex coordinate distributions
1213     TF1* myGaussFit = new TF1("myGaussFit",
1214                               "[0]*exp(-(x-[1])*(x-[1])/(2*[2]*[2]))",
1215                               Vx_Z_Fit->getTH1()->GetXaxis()->GetXmin(),
1216                               Vx_Z_Fit->getTH1()->GetXaxis()->GetXmax());
1217     myGaussFit->SetLineColor(2);
1218     myGaussFit->SetLineWidth(2);
1219     myGaussFit->SetParName(0, "Ampl.");
1220     myGaussFit->SetParName(1, "#mu");
1221     myGaussFit->SetParName(2, "#sigma");
1222 
1223     myGaussFit->SetParameter(0, Vx_X_Fit->getTH1()->GetMaximum());
1224     myGaussFit->SetParameter(1, Vx_X_Fit->getTH1()->GetMean());
1225     myGaussFit->SetParameter(2, Vx_X_Fit->getTH1()->GetRMS());
1226     minXfit = Vx_X_Fit->getTH1()->GetBinLowEdge(1);
1227     for (int i = 0; i < Vx_X_Fit->getTH1()->GetNbinsX(); i++) {
1228       if (Vx_X_Fit->getTH1()->GetBinContent(i + 1) > 0) {
1229         minXfit = Vx_X_Fit->getTH1()->GetBinLowEdge(i + 1);
1230         break;
1231       }
1232     }
1233     maxXfit = Vx_X_Fit->getTH1()->GetBinLowEdge(Vx_X_Fit->getTH1()->GetNbinsX());
1234     for (int i = Vx_X_Fit->getTH1()->GetNbinsX(); i > 0; i--) {
1235       if (Vx_X_Fit->getTH1()->GetBinContent(i) > 0) {
1236         maxXfit = Vx_X_Fit->getTH1()->GetBinLowEdge(i);
1237         break;
1238       }
1239     }
1240     myGaussFit->SetRange(minXfit - (maxXfit - minXfit) / 2., maxXfit + (maxXfit - minXfit) / 2.);
1241     if (Vx_X_Fit->getTH1()->GetEntries() > 0)
1242       Vx_X_Fit->getTH1()->Fit(myGaussFit, "QRL");
1243 
1244     myGaussFit->SetParameter(0, Vx_Y_Fit->getTH1()->GetMaximum());
1245     myGaussFit->SetParameter(1, Vx_Y_Fit->getTH1()->GetMean());
1246     myGaussFit->SetParameter(2, Vx_Y_Fit->getTH1()->GetRMS());
1247     minXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(1);
1248     for (int i = 0; i < Vx_Y_Fit->getTH1()->GetNbinsX(); i++) {
1249       if (Vx_Y_Fit->getTH1()->GetBinContent(i + 1) > 0) {
1250         minXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(i + 1);
1251         break;
1252       }
1253     }
1254     maxXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(Vx_Y_Fit->getTH1()->GetNbinsX());
1255     for (int i = Vx_Y_Fit->getTH1()->GetNbinsX(); i > 0; i--) {
1256       if (Vx_Y_Fit->getTH1()->GetBinContent(i) > 0) {
1257         maxXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(i);
1258         break;
1259       }
1260     }
1261     myGaussFit->SetRange(minXfit - (maxXfit - minXfit) / 2., maxXfit + (maxXfit - minXfit) / 2.);
1262     if (Vx_Y_Fit->getTH1()->GetEntries() > 0)
1263       Vx_Y_Fit->getTH1()->Fit(myGaussFit, "QRL");
1264 
1265     myGaussFit->SetParameter(0, Vx_Z_Fit->getTH1()->GetMaximum());
1266     myGaussFit->SetParameter(1, Vx_Z_Fit->getTH1()->GetMean());
1267     myGaussFit->SetParameter(2, Vx_Z_Fit->getTH1()->GetRMS());
1268     minXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(1);
1269     for (int i = 0; i < Vx_Z_Fit->getTH1()->GetNbinsX(); i++) {
1270       if (Vx_Z_Fit->getTH1()->GetBinContent(i + 1) > 0) {
1271         minXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(i + 1);
1272         break;
1273       }
1274     }
1275     maxXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(Vx_Z_Fit->getTH1()->GetNbinsX());
1276     for (int i = Vx_Z_Fit->getTH1()->GetNbinsX(); i > 0; i--) {
1277       if (Vx_Z_Fit->getTH1()->GetBinContent(i) > 0) {
1278         maxXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(i);
1279         break;
1280       }
1281     }
1282     myGaussFit->SetRange(minXfit - (maxXfit - minXfit) / 2., maxXfit + (maxXfit - minXfit) / 2.);
1283     if (Vx_Z_Fit->getTH1()->GetEntries() > 0)
1284       Vx_Z_Fit->getTH1()->Fit(myGaussFit, "QRL");
1285 
1286     delete myGaussFit;
1287   } else if ((nLumiFit != 0) && (lumiCounter % nLumiFit != 0) && (beginTimeOfFit != 0) && (runNumber != 0)) {
1288     histTitle << "Ongoing: accumulating evts (" << lumiCounter % nLumiFit << " - " << nLumiFit << " in " << lumiCounter
1289               << " - " << maxLumiIntegration << " lumis)";
1290     fitResults->setAxisTitle(histTitle.str(), 1);
1291     if ((debugMode == true) && (outputDebugFile.is_open() == true)) {
1292       outputDebugFile << "\n"
1293                       << "Runnumber " << runNumber << endl;
1294       outputDebugFile << "BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) << " " << (beginTimeOfFit >> 32) << endl;
1295       outputDebugFile << "BeginLumiRange " << beginLumiOfFit << endl;
1296       outputDebugFile << histTitle.str().c_str() << "\n" << endl;
1297     }
1298   } else if ((nLumiFit == 0) || (beginTimeOfFit == 0) || (runNumber == 0)) {
1299     histTitle << "Ongoing: no ongoing fits";
1300     fitResults->setAxisTitle(histTitle.str(), 1);
1301     if ((debugMode == true) && (outputDebugFile.is_open() == true))
1302       outputDebugFile << histTitle.str().c_str() << "\n" << endl;
1303 
1304     endLumiOfFit = lumiBlock.luminosityBlock();
1305 
1306     hitCounter->getTH1()->SetBinContent(endLumiOfFit, (double)totalHits);
1307     hitCounter->getTH1()->SetBinError(endLumiOfFit, std::sqrt((double)totalHits));
1308 
1309     reset("whole");
1310   }
1311 
1312   if (internalDebug == true)
1313     cout << "[Vx3DHLTAnalyzer]::\tHistogram title: " << histTitle.str() << endl;
1314 }
1315 
1316 void Vx3DHLTAnalyzer::bookHistograms(DQMStore::IBooker& ibooker, Run const& iRun, EventSetup const& /* iSetup */) {
1317   ibooker.setCurrentFolder("BeamPixel");
1318 
1319   Vx_X = ibooker.book1D(
1320       "F - vertex x", "Primary Vertex X Distribution", int(rint(xRange / xStep)), -xRange / 2., xRange / 2.);
1321   Vx_Y = ibooker.book1D(
1322       "F - vertex y", "Primary Vertex Y Distribution", int(rint(yRange / yStep)), -yRange / 2., yRange / 2.);
1323   Vx_Z = ibooker.book1D(
1324       "F - vertex z", "Primary Vertex Z Distribution", int(rint(zRange / zStep)), -zRange / 2., zRange / 2.);
1325   Vx_X->setAxisTitle("Primary Vertices X [cm]", 1);
1326   Vx_X->setAxisTitle("Entries [#]", 2);
1327   Vx_Y->setAxisTitle("Primary Vertices Y [cm]", 1);
1328   Vx_Y->setAxisTitle("Entries [#]", 2);
1329   Vx_Z->setAxisTitle("Primary Vertices Z [cm]", 1);
1330   Vx_Z->setAxisTitle("Entries [#]", 2);
1331 
1332   Vx_X_Fit = ibooker.book1D("G - vertex x fit",
1333                             "Primary Vertex X Distribution (For Fit)",
1334                             int(rint(xRange / xStep)),
1335                             -xRange / 2.,
1336                             xRange / 2.);
1337   Vx_Y_Fit = ibooker.book1D("G - vertex y fit",
1338                             "Primary Vertex Y Distribution (For Fit)",
1339                             int(rint(yRange / yStep)),
1340                             -yRange / 2.,
1341                             yRange / 2.);
1342   Vx_Z_Fit = ibooker.book1D("G - vertex z fit",
1343                             "Primary Vertex Z Distribution (For Fit)",
1344                             int(rint(zRange / zStep)),
1345                             -zRange / 2.,
1346                             zRange / 2.);
1347   Vx_X_Fit->setAxisTitle("Primary Vertices X [cm]", 1);
1348   Vx_X_Fit->setAxisTitle("Entries [#]", 2);
1349   Vx_Y_Fit->setAxisTitle("Primary Vertices Y [cm]", 1);
1350   Vx_Y_Fit->setAxisTitle("Entries [#]", 2);
1351   Vx_Z_Fit->setAxisTitle("Primary Vertices Z [cm]", 1);
1352   Vx_Z_Fit->setAxisTitle("Entries [#]", 2);
1353 
1354   Vx_X_Cum = ibooker.book1D("I - vertex x cum",
1355                             "Primary Vertex X Distribution (Cumulative)",
1356                             int(rint(xRange / xStep)),
1357                             -xRange / 2.,
1358                             xRange / 2.);
1359   Vx_Y_Cum = ibooker.book1D("I - vertex y cum",
1360                             "Primary Vertex Y Distribution (Cumulative)",
1361                             int(rint(yRange / yStep)),
1362                             -yRange / 2.,
1363                             yRange / 2.);
1364   Vx_Z_Cum = ibooker.book1D("I - vertex z cum",
1365                             "Primary Vertex Z Distribution (Cumulative)",
1366                             int(rint(zRange / zStep)),
1367                             -zRange / 2.,
1368                             zRange / 2.);
1369   Vx_X_Cum->setAxisTitle("Primary Vertices X [cm]", 1);
1370   Vx_X_Cum->setAxisTitle("Entries [#]", 2);
1371   Vx_Y_Cum->setAxisTitle("Primary Vertices Y [cm]", 1);
1372   Vx_Y_Cum->setAxisTitle("Entries [#]", 2);
1373   Vx_Z_Cum->setAxisTitle("Primary Vertices Z [cm]", 1);
1374   Vx_Z_Cum->setAxisTitle("Entries [#]", 2);
1375 
1376   mXlumi = ibooker.book1D(
1377       "B - muX vs lumi", "#mu_{x} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1378   mYlumi = ibooker.book1D(
1379       "B - muY vs lumi", "#mu_{y} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1380   mZlumi = ibooker.book1D(
1381       "B - muZ vs lumi", "#mu_{z} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1382   mXlumi->setAxisTitle("Lumisection [#]", 1);
1383   mXlumi->setAxisTitle("#mu_{x} [cm]", 2);
1384   mXlumi->getTH1()->SetOption("E1");
1385   mYlumi->setAxisTitle("Lumisection [#]", 1);
1386   mYlumi->setAxisTitle("#mu_{y} [cm]", 2);
1387   mYlumi->getTH1()->SetOption("E1");
1388   mZlumi->setAxisTitle("Lumisection [#]", 1);
1389   mZlumi->setAxisTitle("#mu_{z} [cm]", 2);
1390   mZlumi->getTH1()->SetOption("E1");
1391 
1392   sXlumi = ibooker.book1D(
1393       "C - sigmaX vs lumi", "#sigma_{x} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1394   sYlumi = ibooker.book1D(
1395       "C - sigmaY vs lumi", "#sigma_{y} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1396   sZlumi = ibooker.book1D(
1397       "C - sigmaZ vs lumi", "#sigma_{z} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1398   sXlumi->setAxisTitle("Lumisection [#]", 1);
1399   sXlumi->setAxisTitle("#sigma_{x} [cm]", 2);
1400   sXlumi->getTH1()->SetOption("E1");
1401   sYlumi->setAxisTitle("Lumisection [#]", 1);
1402   sYlumi->setAxisTitle("#sigma_{y} [cm]", 2);
1403   sYlumi->getTH1()->SetOption("E1");
1404   sZlumi->setAxisTitle("Lumisection [#]", 1);
1405   sZlumi->setAxisTitle("#sigma_{z} [cm]", 2);
1406   sZlumi->getTH1()->SetOption("E1");
1407 
1408   dxdzlumi = ibooker.book1D(
1409       "D - dxdz vs lumi", "dX/dZ vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1410   dydzlumi = ibooker.book1D(
1411       "D - dydz vs lumi", "dY/dZ vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1412   dxdzlumi->setAxisTitle("Lumisection [#]", 1);
1413   dxdzlumi->setAxisTitle("dX/dZ [rad]", 2);
1414   dxdzlumi->getTH1()->SetOption("E1");
1415   dydzlumi->setAxisTitle("Lumisection [#]", 1);
1416   dydzlumi->setAxisTitle("dY/dZ [rad]", 2);
1417   dydzlumi->getTH1()->SetOption("E1");
1418 
1419   Vx_ZX = ibooker.book2D("E - vertex zx",
1420                          "Primary Vertex ZX Distribution",
1421                          int(rint(zRange / zStep)),
1422                          -zRange / 2.,
1423                          zRange / 2.,
1424                          int(rint(xRange / xStep)),
1425                          -xRange / 2.,
1426                          xRange / 2.);
1427   Vx_ZY = ibooker.book2D("E - vertex zy",
1428                          "Primary Vertex ZY Distribution",
1429                          int(rint(zRange / zStep)),
1430                          -zRange / 2.,
1431                          zRange / 2.,
1432                          int(rint(yRange / yStep)),
1433                          -yRange / 2.,
1434                          yRange / 2.);
1435   Vx_XY = ibooker.book2D("E - vertex xy",
1436                          "Primary Vertex XY Distribution",
1437                          int(rint(xRange / xStep)),
1438                          -xRange / 2.,
1439                          xRange / 2.,
1440                          int(rint(yRange / yStep)),
1441                          -yRange / 2.,
1442                          yRange / 2.);
1443   Vx_ZX->setAxisTitle("Primary Vertices Z [cm]", 1);
1444   Vx_ZX->setAxisTitle("Primary Vertices X [cm]", 2);
1445   Vx_ZX->setAxisTitle("Entries [#]", 3);
1446   Vx_ZY->setAxisTitle("Primary Vertices Z [cm]", 1);
1447   Vx_ZY->setAxisTitle("Primary Vertices Y [cm]", 2);
1448   Vx_ZY->setAxisTitle("Entries [#]", 3);
1449   Vx_XY->setAxisTitle("Primary Vertices X [cm]", 1);
1450   Vx_XY->setAxisTitle("Primary Vertices Y [cm]", 2);
1451   Vx_XY->setAxisTitle("Entries [#]", 3);
1452 
1453   Vx_ZX_Cum = ibooker.book2D("H - vertex zx cum",
1454                              "Primary Vertex ZX Distribution (Cumulative)",
1455                              int(rint(zRange / zStep)),
1456                              -zRange / 2.,
1457                              zRange / 2.,
1458                              int(rint(xRange / xStep)),
1459                              -xRange / 2.,
1460                              xRange / 2.);
1461   Vx_ZY_Cum = ibooker.book2D("H - vertex zy cum",
1462                              "Primary Vertex ZY Distribution (Cumulative)",
1463                              int(rint(zRange / zStep)),
1464                              -zRange / 2.,
1465                              zRange / 2.,
1466                              int(rint(yRange / yStep)),
1467                              -yRange / 2.,
1468                              yRange / 2.);
1469   Vx_XY_Cum = ibooker.book2D("H - vertex xy cum",
1470                              "Primary Vertex XY Distribution (Cumulative)",
1471                              int(rint(xRange / xStep)),
1472                              -xRange / 2.,
1473                              xRange / 2.,
1474                              int(rint(yRange / yStep)),
1475                              -yRange / 2.,
1476                              yRange / 2.);
1477   Vx_ZX_Cum->setAxisTitle("Primary Vertices Z [cm]", 1);
1478   Vx_ZX_Cum->setAxisTitle("Primary Vertices X [cm]", 2);
1479   Vx_ZX_Cum->setAxisTitle("Entries [#]", 3);
1480   Vx_ZY_Cum->setAxisTitle("Primary Vertices Z [cm]", 1);
1481   Vx_ZY_Cum->setAxisTitle("Primary Vertices Y [cm]", 2);
1482   Vx_ZY_Cum->setAxisTitle("Entries [#]", 3);
1483   Vx_XY_Cum->setAxisTitle("Primary Vertices X [cm]", 1);
1484   Vx_XY_Cum->setAxisTitle("Primary Vertices Y [cm]", 2);
1485   Vx_XY_Cum->setAxisTitle("Entries [#]", 3);
1486 
1487   hitCounter = ibooker.book1D(
1488       "J - pixelHits vs lumi", "# Pixel-Hits vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1489   hitCounter->setAxisTitle("Lumisection [#]", 1);
1490   hitCounter->setAxisTitle("Pixel-Hits [#]", 2);
1491   hitCounter->getTH1()->SetOption("E1");
1492 
1493   goodVxCounter = ibooker.book1D("K - good vertices vs lumi",
1494                                  "# Good vertices vs. Lumisection",
1495                                  nLumiXaxisRange,
1496                                  0.5,
1497                                  ((double)nLumiXaxisRange) + 0.5);
1498   goodVxCounter->setAxisTitle("Lumisection [#]", 1);
1499   goodVxCounter->setAxisTitle("Good vertices [#]", 2);
1500   goodVxCounter->getTH1()->SetOption("E1");
1501 
1502   statusCounter = ibooker.book1D(
1503       "L - status vs lumi", "App. Status vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1504   statusCounter->setAxisTitle("Lumisection [#]", 1);
1505   statusCounter->getTH1()->SetOption("E1");
1506   statusCounter->getTH1()->GetYaxis()->Set(11, -5.5, 5.5);
1507   statusCounter->getTH1()->GetYaxis()->SetBinLabel(1, "Max Lumi.");
1508   statusCounter->getTH1()->GetYaxis()->SetBinLabel(2, "Neg. det.");
1509   statusCounter->getTH1()->GetYaxis()->SetBinLabel(3, "Infinite err.");
1510   statusCounter->getTH1()->GetYaxis()->SetBinLabel(4, "No vtx.");
1511   statusCounter->getTH1()->GetYaxis()->SetBinLabel(5, "Infinite EDM");
1512   statusCounter->getTH1()->GetYaxis()->SetBinLabel(6, "OK");
1513   statusCounter->getTH1()->GetYaxis()->SetBinLabel(7, "MINUIT stat.");
1514   statusCounter->getTH1()->GetYaxis()->SetBinLabel(8, "MINUIT stat.");
1515   statusCounter->getTH1()->GetYaxis()->SetBinLabel(9, "MINUIT stat.");
1516   statusCounter->getTH1()->GetYaxis()->SetBinLabel(10, "MINUIT stat.");
1517   statusCounter->getTH1()->GetYaxis()->SetBinLabel(11, "MINUIT stat.");
1518 
1519   fitResults = ibooker.book2D("A - fit results", "Results of Beam Spot Fit", 2, 0., 2., 9, 0., 9.);
1520   fitResults->setAxisTitle("Ongoing: bootstrapping", 1);
1521   fitResults->setBinLabel(9, "X[cm]", 2);
1522   fitResults->setBinLabel(8, "Y[cm]", 2);
1523   fitResults->setBinLabel(7, "Z[cm]", 2);
1524   fitResults->setBinLabel(6, "#sigma_{X}[cm]", 2);
1525   fitResults->setBinLabel(5, "#sigma_{Y}[cm]", 2);
1526   fitResults->setBinLabel(4, "#sigma_{Z}[cm]", 2);
1527   fitResults->setBinLabel(3, "#frac{dX}{dZ}[rad]", 2);
1528   fitResults->setBinLabel(2, "#frac{dY}{dZ}[rad]", 2);
1529   fitResults->setBinLabel(1, "Vtx[#]", 2);
1530   fitResults->setBinLabel(1, "Value", 1);
1531   fitResults->setBinLabel(2, "Error (stat)", 1);
1532   fitResults->getTH1()->SetOption("text");
1533 
1534   ibooker.setCurrentFolder("BeamPixel/EventInfo");
1535 
1536   reportSummary = ibooker.bookFloat("reportSummary");
1537   reportSummary->Fill(-1);
1538   reportSummaryMap = ibooker.book2D("reportSummaryMap", "Pixel-Vertices Beam Spot: % Good Fits", 1, 0., 1., 1, 0., 1.);
1539   reportSummaryMap->getTH1()->SetBinContent(1, 1, -1);
1540 
1541   ibooker.setCurrentFolder("BeamPixel/EventInfo/reportSummaryContents");
1542 
1543   // Convention for reportSummary and reportSummaryMap:
1544   // - -1%  at the moment of creation of the histogram (i.e. white histogram)
1545   // - n%  numberGoodFits / numberFits
1546 
1547   reset("scratch");  // Initialize histograms after creation
1548 }
1549 
1550 // Define this as a plug-in
1551 DEFINE_FWK_MODULE(Vx3DHLTAnalyzer);