Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:06:53 

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       edm::LogInfo("Vx3DHLTAnalyzer") << "\tI found " << totalHits << " pixel hits until now";
0108       edm::LogInfo("Vx3DHLTAnalyzer") << "\tIn this event there are " << Vx3DCollection->size() << " vertex cadidates";
0109     }
0110 
0111     for (vector<Vertex>::const_iterator it3DVx = Vx3DCollection->begin(); it3DVx != Vx3DCollection->end(); it3DVx++) {
0112       if (internalDebug == true) {
0113         edm::LogInfo("Vx3DHLTAnalyzer") << "\tVertex selections:";
0114         edm::LogInfo("Vx3DHLTAnalyzer") << "\tEvent ID = " << iEvent.id();
0115         edm::LogInfo("Vx3DHLTAnalyzer") << "\tVertex number = " << it3DVx - Vx3DCollection->begin();
0116         edm::LogInfo("Vx3DHLTAnalyzer") << "\tisValid = " << it3DVx->isValid();
0117         edm::LogInfo("Vx3DHLTAnalyzer") << "\tisFake = " << it3DVx->isFake();
0118         edm::LogInfo("Vx3DHLTAnalyzer") << "\tnodof = " << it3DVx->ndof();
0119         edm::LogInfo("Vx3DHLTAnalyzer") << "\ttracksSize = " << it3DVx->tracksSize();
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             edm::LogInfo("Vx3DHLTAnalyzer") << "\tVertex accepted !";
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           edm::LogInfo("Vx3DHLTAnalyzer") << "\tVertex discarded !";
0170 
0171           for (i = 0; i < DIM; i++)
0172             for (j = 0; j < DIM; j++)
0173               edm::LogInfo("Vx3DHLTAnalyzer") << "(i,j) --> " << i << "," << j << " --> " << MyVertex.Covariance[i][j];
0174         }
0175       } else if (internalDebug == true)
0176         edm::LogInfo("Vx3DHLTAnalyzer") << "\tVertex discarded !";
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 == NOT OK - fit status (MINUIT manual)   #
0271   // #  0 == OK                                    #
0272   // # -1 == NOT OK - not finite edm               #
0273   // # -2 == NOT OK - not enough "minNentries"     #
0274   // # -3 == NOT OK - not finite errors            #
0275   // # -4 == NOT OK - negative determinant         #
0276   // # -5 == NOT OK - maxLumiIntegration reached   #
0277   // # -6 == NOT OK - FatalRootError: @SUB=Minuit2 #
0278   // ###############################################
0279 
0280   if ((vals != nullptr) && (vals->size() == nParams * 2)) {
0281     double nSigmaXY = 10.;
0282     double nSigmaZ = 10.;
0283     double parDistanceXY = 1e-3;   // Unit: [cm]
0284     double parDistanceZ = 1e-2;    // Unit: [cm]
0285     double parDistanceddZ = 1e-3;  // Unit: [rad]
0286     double parDistanceCxy = 1e-5;  // Unit: [cm^2]
0287     double bestEdm;
0288 
0289     const unsigned int trials = 4;
0290     double largerDist[trials] = {0.1, 5., 10., 100.};
0291 
0292     double covxz, covyz, det;
0293     double deltaMean;
0294     int bestMovementX = 1;
0295     int bestMovementY = 1;
0296     int bestMovementZ = 1;
0297     int goodData;
0298 
0299     double edm;
0300 
0301     vector<double>::const_iterator it = vals->begin();
0302 
0303     ROOT::Math::Minimizer* Gauss3D = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
0304     Gauss3D->SetErrorDef(1.0);
0305     if (internalDebug == true)
0306       Gauss3D->SetPrintLevel(3);
0307     else
0308       Gauss3D->SetPrintLevel(0);
0309 
0310     ROOT::Math::Functor _Gauss3DFunc(this, &Vx3DHLTAnalyzer::Gauss3DFunc, nParams);
0311     Gauss3D->SetFunction(_Gauss3DFunc);
0312 
0313     if (internalDebug == true)
0314       edm::LogInfo("Vx3DHLTAnalyzer") << "\t@@@ START FITTING @@@";
0315 
0316     // @@@ Fit at X-deltaMean | X | X+deltaMean @@@
0317     bestEdm = 1.;
0318     for (int i = 0; i < 3; i++) {
0319       deltaMean = (double(i) - 1.) * std::sqrt(*(it + 0));
0320       if (internalDebug == true)
0321         edm::LogInfo("Vx3DHLTAnalyzer") << "\tdeltaMean --> " << deltaMean;
0322 
0323       Gauss3D->Clear();
0324 
0325       Gauss3D->SetVariable(0, "var x ", *(it + 0), parDistanceXY * parDistanceXY);
0326       Gauss3D->SetVariable(1, "var y ", *(it + 1), parDistanceXY * parDistanceXY);
0327       Gauss3D->SetVariable(2, "var z ", *(it + 2), parDistanceZ * parDistanceZ);
0328       Gauss3D->SetVariable(3, "cov xy", *(it + 3), parDistanceCxy);
0329       Gauss3D->SetVariable(4, "dydz  ", *(it + 4), parDistanceddZ);
0330       Gauss3D->SetVariable(5, "dxdz  ", *(it + 5), parDistanceddZ);
0331       Gauss3D->SetVariable(6, "mean x", *(it + 6) + deltaMean, parDistanceXY);
0332       Gauss3D->SetVariable(7, "mean y", *(it + 7), parDistanceXY);
0333       Gauss3D->SetVariable(8, "mean z", *(it + 8), parDistanceZ);
0334 
0335       // Set the central positions of the centroid for vertex rejection
0336       xPos = *(it + 6) + deltaMean;
0337       yPos = *(it + 7);
0338       zPos = *(it + 8);
0339 
0340       // Set dimensions of the centroid for vertex rejection
0341       maxTransRadius = nSigmaXY * std::sqrt(std::fabs(*(it + 0)) + std::fabs(*(it + 1)) / 2.);
0342       maxLongLength = nSigmaZ * std::sqrt(std::fabs(*(it + 2)));
0343 
0344       try {
0345         Gauss3D->Minimize();
0346         goodData = Gauss3D->Status();
0347         edm = Gauss3D->Edm();
0348       } catch (cms::Exception& er) {
0349         edm::LogError("Vx3DHLTAnalyzer") << "\tCaught Minuit2 exception @ Fit at X: \n" << er.what();
0350         goodData = -6;
0351         edm = 1.;
0352         continue;
0353       }
0354 
0355       if (counterVx < minNentries)
0356         goodData = -2;
0357       else if (isNotFinite(edm) == true) {
0358         goodData = -1;
0359         if (internalDebug == true)
0360           edm::LogInfo("Vx3DHLTAnalyzer") << "\tNot finite edm !";
0361       } else if (goodData != -6)
0362         for (unsigned int j = 0; j < nParams; j++)
0363           if (isNotFinite(Gauss3D->Errors()[j]) == true) {
0364             goodData = -3;
0365             if (internalDebug == true)
0366               edm::LogInfo("Vx3DHLTAnalyzer") << "\tNot finite errors !";
0367             break;
0368           }
0369       if (goodData == 0) {
0370         covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
0371                 Gauss3D->X()[5] * Gauss3D->X()[3];
0372         covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
0373                 Gauss3D->X()[4] * Gauss3D->X()[3];
0374 
0375         det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
0376               Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
0377               covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
0378         if (det < 0.) {
0379           goodData = -4;
0380           if (internalDebug == true)
0381             edm::LogInfo("Vx3DHLTAnalyzer") << "\tNegative determinant !";
0382         }
0383       }
0384 
0385       if ((goodData == 0) && (std::fabs(edm) < bestEdm)) {
0386         bestEdm = edm;
0387         bestMovementX = i;
0388       }
0389     }
0390     if (internalDebug == true)
0391       edm::LogInfo("Vx3DHLTAnalyzer") << "\tFound bestMovementX --> " << bestMovementX;
0392 
0393     // @@@ Fit at Y-deltaMean | Y | Y+deltaMean @@@
0394     bestEdm = 1.;
0395     for (int i = 0; i < 3; i++) {
0396       deltaMean = (double(i) - 1.) * std::sqrt(*(it + 1));
0397       if (internalDebug == true) {
0398         edm::LogInfo("Vx3DHLTAnalyzer") << "\tdeltaMean --> " << deltaMean;
0399         edm::LogInfo("Vx3DHLTAnalyzer") << "\tdeltaMean X --> " << (double(bestMovementX) - 1.) * std::sqrt(*(it + 0));
0400       }
0401 
0402       Gauss3D->Clear();
0403 
0404       Gauss3D->SetVariable(0, "var x ", *(it + 0), parDistanceXY * parDistanceXY);
0405       Gauss3D->SetVariable(1, "var y ", *(it + 1), parDistanceXY * parDistanceXY);
0406       Gauss3D->SetVariable(2, "var z ", *(it + 2), parDistanceZ * parDistanceZ);
0407       Gauss3D->SetVariable(3, "cov xy", *(it + 3), parDistanceCxy);
0408       Gauss3D->SetVariable(4, "dydz  ", *(it + 4), parDistanceddZ);
0409       Gauss3D->SetVariable(5, "dxdz  ", *(it + 5), parDistanceddZ);
0410       Gauss3D->SetVariable(6, "mean x", *(it + 6) + (double(bestMovementX) - 1.) * std::sqrt(*(it + 0)), parDistanceXY);
0411       Gauss3D->SetVariable(7, "mean y", *(it + 7) + deltaMean, parDistanceXY);
0412       Gauss3D->SetVariable(8, "mean z", *(it + 8), parDistanceZ);
0413 
0414       // Set the central positions of the centroid for vertex rejection
0415       xPos = *(it + 6) + (double(bestMovementX) - 1.) * std::sqrt(*(it + 0));
0416       yPos = *(it + 7) + deltaMean;
0417       zPos = *(it + 8);
0418 
0419       // Set dimensions of the centroid for vertex rejection
0420       maxTransRadius = nSigmaXY * std::sqrt(std::fabs(*(it + 0)) + std::fabs(*(it + 1)) / 2.);
0421       maxLongLength = nSigmaZ * std::sqrt(std::fabs(*(it + 2)));
0422 
0423       try {
0424         Gauss3D->Minimize();
0425         goodData = Gauss3D->Status();
0426         edm = Gauss3D->Edm();
0427       } catch (cms::Exception& er) {
0428         edm::LogError("Vx3DHLTAnalyzer") << "\tCaught Minuit2 exception @ Fit at Y: \n" << er.what();
0429         goodData = -6;
0430         edm = 1.;
0431         continue;
0432       }
0433 
0434       if (counterVx < minNentries)
0435         goodData = -2;
0436       else if (isNotFinite(edm) == true) {
0437         goodData = -1;
0438         if (internalDebug == true)
0439           edm::LogInfo("Vx3DHLTAnalyzer") << "\tNot finite edm !";
0440       } else if (goodData != -6)
0441         for (unsigned int j = 0; j < nParams; j++)
0442           if (isNotFinite(Gauss3D->Errors()[j]) == true) {
0443             goodData = -3;
0444             if (internalDebug == true)
0445               edm::LogInfo("Vx3DHLTAnalyzer") << "\tNot finite errors !";
0446             break;
0447           }
0448       if (goodData == 0) {
0449         covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
0450                 Gauss3D->X()[5] * Gauss3D->X()[3];
0451         covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
0452                 Gauss3D->X()[4] * Gauss3D->X()[3];
0453 
0454         det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
0455               Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
0456               covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
0457         if (det < 0.) {
0458           goodData = -4;
0459           if (internalDebug == true)
0460             edm::LogInfo("Vx3DHLTAnalyzer") << "\tNegative determinant !";
0461         }
0462       }
0463 
0464       if ((goodData == 0) && (std::fabs(edm) < bestEdm)) {
0465         bestEdm = edm;
0466         bestMovementY = i;
0467       }
0468     }
0469     if (internalDebug == true)
0470       edm::LogInfo("Vx3DHLTAnalyzer") << "\tFound bestMovementY --> " << bestMovementY;
0471 
0472     // @@@ Fit at Z-deltaMean | Z | Z+deltaMean @@@
0473     bestEdm = 1.;
0474     for (int i = 0; i < 3; i++) {
0475       deltaMean = (double(i) - 1.) * std::sqrt(*(it + 2));
0476       if (internalDebug == true) {
0477         edm::LogInfo("Vx3DHLTAnalyzer") << "\tdeltaMean --> " << deltaMean;
0478         edm::LogInfo("Vx3DHLTAnalyzer") << "\tdeltaMean X --> " << (double(bestMovementX) - 1.) * std::sqrt(*(it + 0));
0479         edm::LogInfo("Vx3DHLTAnalyzer") << "\tdeltaMean Y --> " << (double(bestMovementY) - 1.) * std::sqrt(*(it + 1));
0480       }
0481 
0482       Gauss3D->Clear();
0483 
0484       Gauss3D->SetVariable(0, "var x ", *(it + 0), parDistanceXY * parDistanceXY);
0485       Gauss3D->SetVariable(1, "var y ", *(it + 1), parDistanceXY * parDistanceXY);
0486       Gauss3D->SetVariable(2, "var z ", *(it + 2), parDistanceZ * parDistanceZ);
0487       Gauss3D->SetVariable(3, "cov xy", *(it + 3), parDistanceCxy);
0488       Gauss3D->SetVariable(4, "dydz  ", *(it + 4), parDistanceddZ);
0489       Gauss3D->SetVariable(5, "dxdz  ", *(it + 5), parDistanceddZ);
0490       Gauss3D->SetVariable(6, "mean x", *(it + 6) + (double(bestMovementX) - 1.) * std::sqrt(*(it + 0)), parDistanceXY);
0491       Gauss3D->SetVariable(7, "mean y", *(it + 7) + (double(bestMovementY) - 1.) * std::sqrt(*(it + 1)), parDistanceXY);
0492       Gauss3D->SetVariable(8, "mean z", *(it + 8) + deltaMean, parDistanceZ);
0493 
0494       // Set the central positions of the centroid for vertex rejection
0495       xPos = *(it + 6) + (double(bestMovementX) - 1.) * std::sqrt(*(it + 0));
0496       yPos = *(it + 7) + (double(bestMovementY) - 1.) * std::sqrt(*(it + 1));
0497       zPos = *(it + 8) + deltaMean;
0498 
0499       // Set dimensions of the centroid for vertex rejection
0500       maxTransRadius = nSigmaXY * std::sqrt(std::fabs(*(it + 0)) + std::fabs(*(it + 1)) / 2.);
0501       maxLongLength = nSigmaZ * std::sqrt(std::fabs(*(it + 2)));
0502 
0503       try {
0504         Gauss3D->Minimize();
0505         goodData = Gauss3D->Status();
0506         edm = Gauss3D->Edm();
0507       } catch (cms::Exception& er) {
0508         edm::LogError("Vx3DHLTAnalyzer") << "\tCaught Minuit2 exception @ Fit at Z: \n" << er.what();
0509         goodData = -6;
0510         edm = 1.;
0511         continue;
0512       }
0513 
0514       if (counterVx < minNentries)
0515         goodData = -2;
0516       else if (isNotFinite(edm) == true) {
0517         goodData = -1;
0518         if (internalDebug == true)
0519           edm::LogInfo("Vx3DHLTAnalyzer") << "\tNot finite edm !";
0520       } else if (goodData != -6)
0521         for (unsigned int j = 0; j < nParams; j++)
0522           if (isNotFinite(Gauss3D->Errors()[j]) == true) {
0523             goodData = -3;
0524             if (internalDebug == true)
0525               edm::LogInfo("Vx3DHLTAnalyzer") << "\tNot finite errors !";
0526             break;
0527           }
0528       if (goodData == 0) {
0529         covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
0530                 Gauss3D->X()[5] * Gauss3D->X()[3];
0531         covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
0532                 Gauss3D->X()[4] * Gauss3D->X()[3];
0533 
0534         det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
0535               Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
0536               covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
0537         if (det < 0.) {
0538           goodData = -4;
0539           if (internalDebug == true)
0540             edm::LogInfo("Vx3DHLTAnalyzer") << "\tNegative determinant !";
0541         }
0542       }
0543 
0544       if ((goodData == 0) && (std::fabs(edm) < bestEdm)) {
0545         bestEdm = edm;
0546         bestMovementZ = i;
0547       }
0548     }
0549     if (internalDebug == true)
0550       edm::LogInfo("Vx3DHLTAnalyzer") << "\tFound bestMovementZ --> " << bestMovementZ;
0551 
0552     Gauss3D->Clear();
0553 
0554     // @@@ FINAL FIT @@@
0555     Gauss3D->SetVariable(0, "var x ", *(it + 0), parDistanceXY * parDistanceXY);
0556     Gauss3D->SetVariable(1, "var y ", *(it + 1), parDistanceXY * parDistanceXY);
0557     Gauss3D->SetVariable(2, "var z ", *(it + 2), parDistanceZ * parDistanceZ);
0558     Gauss3D->SetVariable(3, "cov xy", *(it + 3), parDistanceCxy);
0559     Gauss3D->SetVariable(4, "dydz  ", *(it + 4), parDistanceddZ);
0560     Gauss3D->SetVariable(5, "dxdz  ", *(it + 5), parDistanceddZ);
0561     Gauss3D->SetVariable(6, "mean x", *(it + 6) + (double(bestMovementX) - 1.) * std::sqrt(*(it + 0)), parDistanceXY);
0562     Gauss3D->SetVariable(7, "mean y", *(it + 7) + (double(bestMovementY) - 1.) * std::sqrt(*(it + 1)), parDistanceXY);
0563     Gauss3D->SetVariable(8, "mean z", *(it + 8) + (double(bestMovementZ) - 1.) * std::sqrt(*(it + 2)), parDistanceZ);
0564 
0565     // Set the central positions of the centroid for vertex rejection
0566     xPos = *(it + 6) + (double(bestMovementX) - 1.) * std::sqrt(*(it + 0));
0567     yPos = *(it + 7) + (double(bestMovementY) - 1.) * std::sqrt(*(it + 1));
0568     zPos = *(it + 8) + (double(bestMovementZ) - 1.) * std::sqrt(*(it + 2));
0569 
0570     // Set dimensions of the centroid for vertex rejection
0571     maxTransRadius = nSigmaXY * std::sqrt(std::fabs(*(it + 0)) + std::fabs(*(it + 1)) / 2.);
0572     maxLongLength = nSigmaZ * std::sqrt(std::fabs(*(it + 2)));
0573 
0574     try {
0575       Gauss3D->Minimize();
0576       goodData = Gauss3D->Status();
0577       edm = Gauss3D->Edm();
0578     } catch (cms::Exception& er) {
0579       edm::LogError("Vx3DHLTAnalyzer") << "\tCaught Minuit2 exception @ Final fit: \n" << er.what();
0580       goodData = -6;
0581       edm = 1.;
0582     }
0583 
0584     if (counterVx < minNentries)
0585       goodData = -2;
0586     else if (isNotFinite(edm) == true) {
0587       goodData = -1;
0588       if (internalDebug == true)
0589         edm::LogInfo("Vx3DHLTAnalyzer") << "\tNot finite edm !";
0590     } else if (goodData != -6)
0591       for (unsigned int j = 0; j < nParams; j++)
0592         if (isNotFinite(Gauss3D->Errors()[j]) == true) {
0593           goodData = -3;
0594           if (internalDebug == true)
0595             edm::LogInfo("Vx3DHLTAnalyzer") << "\tNot finite errors !";
0596           break;
0597         }
0598     if (goodData == 0) {
0599       covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
0600               Gauss3D->X()[5] * Gauss3D->X()[3];
0601       covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
0602               Gauss3D->X()[4] * Gauss3D->X()[3];
0603 
0604       det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
0605             Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
0606             covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
0607       if (det < 0.) {
0608         goodData = -4;
0609         if (internalDebug == true)
0610           edm::LogInfo("Vx3DHLTAnalyzer") << "\tNegative determinant !";
0611       }
0612     }
0613 
0614     // @@@ FIT WITH DIFFERENT PARAMETER DISTANCES @@@
0615     for (unsigned int i = 0; i < trials; i++) {
0616       if ((goodData != 0) && (goodData != -2)) {
0617         Gauss3D->Clear();
0618 
0619         if (internalDebug == true)
0620           edm::LogInfo("Vx3DHLTAnalyzer") << "\tFIT WITH DIFFERENT PARAMETER DISTANCES - STEP " << i + 1;
0621 
0622         Gauss3D->SetVariable(0, "var x ", *(it + 0), parDistanceXY * parDistanceXY * largerDist[i]);
0623         Gauss3D->SetVariable(1, "var y ", *(it + 1), parDistanceXY * parDistanceXY * largerDist[i]);
0624         Gauss3D->SetVariable(2, "var z ", *(it + 2), parDistanceZ * parDistanceZ * largerDist[i]);
0625         Gauss3D->SetVariable(3, "cov xy", *(it + 3), parDistanceCxy * largerDist[i]);
0626         Gauss3D->SetVariable(4, "dydz  ", *(it + 4), parDistanceddZ * largerDist[i]);
0627         Gauss3D->SetVariable(5, "dxdz  ", *(it + 5), parDistanceddZ * largerDist[i]);
0628         Gauss3D->SetVariable(6,
0629                              "mean x",
0630                              *(it + 6) + (double(bestMovementX) - 1.) * std::sqrt(*(it + 0)),
0631                              parDistanceXY * largerDist[i]);
0632         Gauss3D->SetVariable(7,
0633                              "mean y",
0634                              *(it + 7) + (double(bestMovementY) - 1.) * std::sqrt(*(it + 1)),
0635                              parDistanceXY * largerDist[i]);
0636         Gauss3D->SetVariable(
0637             8, "mean z", *(it + 8) + (double(bestMovementZ) - 1.) * std::sqrt(*(it + 2)), parDistanceZ * largerDist[i]);
0638 
0639         // Set the central positions of the centroid for vertex rejection
0640         xPos = *(it + 6) + (double(bestMovementX) - 1.) * std::sqrt(*(it + 0));
0641         yPos = *(it + 7) + (double(bestMovementY) - 1.) * std::sqrt(*(it + 1));
0642         zPos = *(it + 8) + (double(bestMovementZ) - 1.) * std::sqrt(*(it + 2));
0643 
0644         // Set dimensions of the centroid for vertex rejection
0645         maxTransRadius = nSigmaXY * std::sqrt(std::fabs(*(it + 0)) + std::fabs(*(it + 1)) / 2.);
0646         maxLongLength = nSigmaZ * std::sqrt(std::fabs(*(it + 2)));
0647 
0648         try {
0649           Gauss3D->Minimize();
0650           goodData = Gauss3D->Status();
0651           edm = Gauss3D->Edm();
0652         } catch (cms::Exception& er) {
0653           edm::LogError("Vx3DHLTAnalyzer") << "\tCaught Minuit2 exception @ Fit with different distances: \n"
0654                                            << er.what();
0655           goodData = -6;
0656           edm = 1.;
0657           continue;
0658         }
0659 
0660         if (counterVx < minNentries)
0661           goodData = -2;
0662         else if (isNotFinite(edm) == true) {
0663           goodData = -1;
0664           if (internalDebug == true)
0665             edm::LogInfo("Vx3DHLTAnalyzer") << "\tNot finite edm !";
0666         } else if (goodData != -6)
0667           for (unsigned int j = 0; j < nParams; j++)
0668             if (isNotFinite(Gauss3D->Errors()[j]) == true) {
0669               goodData = -3;
0670               if (internalDebug == true)
0671                 edm::LogInfo("Vx3DHLTAnalyzer") << "\tNot finite errors !";
0672               break;
0673             }
0674         if (goodData == 0) {
0675           covyz = Gauss3D->X()[4] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[1])) -
0676                   Gauss3D->X()[5] * Gauss3D->X()[3];
0677           covxz = Gauss3D->X()[5] * (std::fabs(Gauss3D->X()[2]) - std::fabs(Gauss3D->X()[0])) -
0678                   Gauss3D->X()[4] * Gauss3D->X()[3];
0679 
0680           det = std::fabs(Gauss3D->X()[0]) * (std::fabs(Gauss3D->X()[1]) * std::fabs(Gauss3D->X()[2]) - covyz * covyz) -
0681                 Gauss3D->X()[3] * (Gauss3D->X()[3] * std::fabs(Gauss3D->X()[2]) - covxz * covyz) +
0682                 covxz * (Gauss3D->X()[3] * covyz - covxz * std::fabs(Gauss3D->X()[1]));
0683           if (det < 0.) {
0684             goodData = -4;
0685             if (internalDebug == true)
0686               edm::LogInfo("Vx3DHLTAnalyzer") << "\tNegative determinant !";
0687           }
0688         }
0689       } else
0690         break;
0691     }
0692 
0693     if (goodData == 0)
0694       for (unsigned int i = 0; i < nParams; i++) {
0695         vals->operator[](i) = Gauss3D->X()[i];
0696         vals->operator[](i + nParams) = Gauss3D->Errors()[i];
0697       }
0698 
0699     delete Gauss3D;
0700     return goodData;
0701   }
0702 
0703   return -1;
0704 }
0705 
0706 void Vx3DHLTAnalyzer::reset(string ResetType) {
0707   if ((debugMode == true) && (outputDebugFile.is_open() == true)) {
0708     outputDebugFile << "Runnumber " << runNumber << endl;
0709     outputDebugFile << "BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) << " " << (beginTimeOfFit >> 32) << endl;
0710     outputDebugFile << "BeginLumiRange " << beginLumiOfFit << endl;
0711     outputDebugFile << "EndTimeOfFit " << formatTime(endTimeOfFit >> 32) << " " << (endTimeOfFit >> 32) << endl;
0712     outputDebugFile << "EndLumiRange " << endLumiOfFit << endl;
0713     outputDebugFile << "LumiCounter " << lumiCounter << endl;
0714     outputDebugFile << "LastLumiOfFit " << lastLumiOfFit << endl;
0715   }
0716 
0717   if (ResetType == "scratch") {
0718     runNumber = 0;
0719     numberGoodFits = 0;
0720     numberFits = 0;
0721     lastLumiOfFit = 0;
0722 
0723     Vx_X->Reset();
0724     Vx_Y->Reset();
0725     Vx_Z->Reset();
0726 
0727     Vx_X_Fit->Reset();
0728     Vx_Y_Fit->Reset();
0729     Vx_Z_Fit->Reset();
0730 
0731     Vx_ZX->Reset();
0732     Vx_ZY->Reset();
0733     Vx_XY->Reset();
0734 
0735     Vx_X_Cum->Reset();
0736     Vx_Y_Cum->Reset();
0737     Vx_Z_Cum->Reset();
0738 
0739     Vx_ZX_Cum->Reset();
0740     Vx_ZY_Cum->Reset();
0741     Vx_XY_Cum->Reset();
0742 
0743     mXlumi->Reset();
0744     mYlumi->Reset();
0745     mZlumi->Reset();
0746 
0747     sXlumi->Reset();
0748     sYlumi->Reset();
0749     sZlumi->Reset();
0750 
0751     dxdzlumi->Reset();
0752     dydzlumi->Reset();
0753 
0754     hitCounter->Reset();
0755     goodVxCounter->Reset();
0756     statusCounter->Reset();
0757     fitResults->Reset();
0758 
0759     reportSummary->Fill(-1);
0760     reportSummaryMap->getTH1()->SetBinContent(1, 1, -1);
0761 
0762     Vertices.clear();
0763 
0764     lumiCounter = 0;
0765     totalHits = 0;
0766     beginTimeOfFit = 0;
0767     endTimeOfFit = 0;
0768     beginLumiOfFit = 0;
0769     endLumiOfFit = 0;
0770 
0771     if (internalDebug == true)
0772       edm::LogInfo("Vx3DHLTAnalyzer") << "\tReset issued: scratch";
0773     if ((debugMode == true) && (outputDebugFile.is_open() == true))
0774       outputDebugFile << "Reset -scratch- issued\n" << endl;
0775   } else if (ResetType == "whole") {
0776     Vx_X->Reset();
0777     Vx_Y->Reset();
0778     Vx_Z->Reset();
0779 
0780     Vx_ZX->Reset();
0781     Vx_ZY->Reset();
0782     Vx_XY->Reset();
0783 
0784     Vertices.clear();
0785 
0786     lumiCounter = 0;
0787     totalHits = 0;
0788     beginTimeOfFit = 0;
0789     endTimeOfFit = 0;
0790     beginLumiOfFit = 0;
0791     endLumiOfFit = 0;
0792 
0793     if (internalDebug == true)
0794       edm::LogInfo("Vx3DHLTAnalyzer") << "\tReset issued: whole";
0795     if ((debugMode == true) && (outputDebugFile.is_open() == true))
0796       outputDebugFile << "Reset -whole- issued\n" << endl;
0797   } else if (ResetType == "fit") {
0798     Vx_X_Fit->Reset();
0799     Vx_Y_Fit->Reset();
0800     Vx_Z_Fit->Reset();
0801 
0802     if (internalDebug == true)
0803       edm::LogInfo("Vx3DHLTAnalyzer") << "\tReset issued: fit";
0804     if ((debugMode == true) && (outputDebugFile.is_open() == true))
0805       outputDebugFile << "Reset -fit- issued\n" << endl;
0806   } else if (ResetType == "hitCounter") {
0807     totalHits = 0;
0808 
0809     if (internalDebug == true)
0810       edm::LogInfo("Vx3DHLTAnalyzer") << "\tReset issued: hitCounter";
0811     if ((debugMode == true) && (outputDebugFile.is_open() == true))
0812       outputDebugFile << "Reset -hitCounter- issued\n" << endl;
0813   }
0814 }
0815 
0816 void Vx3DHLTAnalyzer::writeToFile(vector<double>* vals,
0817                                   TimeValue_t BeginTimeOfFit,
0818                                   TimeValue_t EndTimeOfFit,
0819                                   unsigned int BeginLumiOfFit,
0820                                   unsigned int EndLumiOfFit,
0821                                   int dataType) {
0822   stringstream BufferString;
0823   BufferString.precision(5);
0824 
0825   outputFile.open(fileName.c_str(), ios::out);
0826 
0827   if ((outputFile.is_open() == true) && (vals != nullptr) && (vals->size() == (nParams - 1) * 2)) {
0828     vector<double>::const_iterator it = vals->begin();
0829 
0830     outputFile << "Runnumber " << runNumber << endl;
0831     outputFile << "BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) << " " << (beginTimeOfFit >> 32) << endl;
0832     outputFile << "EndTimeOfFit " << formatTime(endTimeOfFit >> 32) << " " << (endTimeOfFit >> 32) << endl;
0833     outputFile << "LumiRange " << beginLumiOfFit << " - " << endLumiOfFit << endl;
0834     outputFile << "Type " << dataType << endl;
0835     // 3D Vertexing with Pixel Tracks:
0836     // Good data = Type  3
0837     // Bad data  = Type -1
0838 
0839     BufferString << *(it + 0);
0840     outputFile << "X0 " << BufferString.str().c_str() << endl;
0841     BufferString.str("");
0842 
0843     BufferString << *(it + 1);
0844     outputFile << "Y0 " << BufferString.str().c_str() << endl;
0845     BufferString.str("");
0846 
0847     BufferString << *(it + 2);
0848     outputFile << "Z0 " << BufferString.str().c_str() << endl;
0849     BufferString.str("");
0850 
0851     BufferString << *(it + 5);
0852     outputFile << "sigmaZ0 " << BufferString.str().c_str() << endl;
0853     BufferString.str("");
0854 
0855     BufferString << *(it + 6);
0856     outputFile << "dxdz " << BufferString.str().c_str() << endl;
0857     BufferString.str("");
0858 
0859     BufferString << *(it + 7);
0860     outputFile << "dydz " << BufferString.str().c_str() << endl;
0861     BufferString.str("");
0862 
0863     BufferString << *(it + 3);
0864     outputFile << "BeamWidthX " << BufferString.str().c_str() << endl;
0865     BufferString.str("");
0866 
0867     BufferString << *(it + 4);
0868     outputFile << "BeamWidthY " << BufferString.str().c_str() << endl;
0869     BufferString.str("");
0870 
0871     outputFile << "Cov(0,j) " << *(it + 8) << " 0 0 0 0 0 0" << endl;
0872     outputFile << "Cov(1,j) 0 " << *(it + 9) << " 0 0 0 0 0" << endl;
0873     outputFile << "Cov(2,j) 0 0 " << *(it + 10) << " 0 0 0 0" << endl;
0874     outputFile << "Cov(3,j) 0 0 0 " << *(it + 13) << " 0 0 0" << endl;
0875     outputFile << "Cov(4,j) 0 0 0 0 " << *(it + 14) << " 0 0" << endl;
0876     outputFile << "Cov(5,j) 0 0 0 0 0 " << *(it + 15) << " 0" << endl;
0877     outputFile << "Cov(6,j) 0 0 0 0 0 0 "
0878                << ((*(it + 11)) + (*(it + 12)) + 2. * std::sqrt((*(it + 11)) * (*(it + 12)))) / 4. << endl;
0879 
0880     outputFile << "EmittanceX 0" << endl;
0881     outputFile << "EmittanceY 0" << endl;
0882     outputFile << "BetaStar 0" << endl;
0883   }
0884   outputFile.close();
0885 
0886   if ((debugMode == true) && (outputDebugFile.is_open() == true) && (vals != nullptr) &&
0887       (vals->size() == (nParams - 1) * 2)) {
0888     vector<double>::const_iterator it = vals->begin();
0889 
0890     outputDebugFile << "Runnumber " << runNumber << endl;
0891     outputDebugFile << "BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) << " " << (beginTimeOfFit >> 32) << endl;
0892     outputDebugFile << "EndTimeOfFit " << formatTime(endTimeOfFit >> 32) << " " << (endTimeOfFit >> 32) << endl;
0893     outputDebugFile << "LumiRange " << beginLumiOfFit << " - " << endLumiOfFit << endl;
0894     outputDebugFile << "Type " << dataType << endl;
0895     // 3D Vertexing with Pixel Tracks:
0896     // Good data = Type  3
0897     // Bad data  = Type -1
0898 
0899     BufferString << *(it + 0);
0900     outputDebugFile << "X0 " << BufferString.str().c_str() << endl;
0901     BufferString.str("");
0902 
0903     BufferString << *(it + 1);
0904     outputDebugFile << "Y0 " << BufferString.str().c_str() << endl;
0905     BufferString.str("");
0906 
0907     BufferString << *(it + 2);
0908     outputDebugFile << "Z0 " << BufferString.str().c_str() << endl;
0909     BufferString.str("");
0910 
0911     BufferString << *(it + 5);
0912     outputDebugFile << "sigmaZ0 " << BufferString.str().c_str() << endl;
0913     BufferString.str("");
0914 
0915     BufferString << *(it + 6);
0916     outputDebugFile << "dxdz " << BufferString.str().c_str() << endl;
0917     BufferString.str("");
0918 
0919     BufferString << *(it + 7);
0920     outputDebugFile << "dydz " << BufferString.str().c_str() << endl;
0921     BufferString.str("");
0922 
0923     BufferString << *(it + 3);
0924     outputDebugFile << "BeamWidthX " << BufferString.str().c_str() << endl;
0925     BufferString.str("");
0926 
0927     BufferString << *(it + 4);
0928     outputDebugFile << "BeamWidthY " << BufferString.str().c_str() << endl;
0929     BufferString.str("");
0930 
0931     outputDebugFile << "Cov(0,j) " << *(it + 8) << " 0 0 0 0 0 0" << endl;
0932     outputDebugFile << "Cov(1,j) 0 " << *(it + 9) << " 0 0 0 0 0" << endl;
0933     outputDebugFile << "Cov(2,j) 0 0 " << *(it + 10) << " 0 0 0 0" << endl;
0934     outputDebugFile << "Cov(3,j) 0 0 0 " << *(it + 13) << " 0 0 0" << endl;
0935     outputDebugFile << "Cov(4,j) 0 0 0 0 " << *(it + 14) << " 0 0" << endl;
0936     outputDebugFile << "Cov(5,j) 0 0 0 0 0 " << *(it + 15) << " 0" << endl;
0937     outputDebugFile << "Cov(6,j) 0 0 0 0 0 0 "
0938                     << ((*(it + 11)) + (*(it + 12)) + 2. * std::sqrt((*(it + 11)) * (*(it + 12)))) / 4. << endl;
0939 
0940     outputDebugFile << "EmittanceX 0" << endl;
0941     outputDebugFile << "EmittanceY 0" << endl;
0942     outputDebugFile << "BetaStar 0" << endl;
0943 
0944     outputDebugFile << "\nUsed vertices: " << counterVx << "\n" << endl;
0945   }
0946 }
0947 
0948 void Vx3DHLTAnalyzer::printFitParams(const vector<double>& fitResults) {
0949   edm::LogInfo("Vx3DHLTAnalyzer") << "var x -->  " << fitResults[0] << " +/- " << fitResults[0 + nParams];
0950   edm::LogInfo("Vx3DHLTAnalyzer") << "var y -->  " << fitResults[1] << " +/- " << fitResults[1 + nParams];
0951   edm::LogInfo("Vx3DHLTAnalyzer") << "var z -->  " << fitResults[2] << " +/- " << fitResults[2 + nParams];
0952   edm::LogInfo("Vx3DHLTAnalyzer") << "cov xy --> " << fitResults[3] << " +/- " << fitResults[3 + nParams];
0953   edm::LogInfo("Vx3DHLTAnalyzer") << "dydz   --> " << fitResults[4] << " +/- " << fitResults[4 + nParams];
0954   edm::LogInfo("Vx3DHLTAnalyzer") << "dxdz   --> " << fitResults[5] << " +/- " << fitResults[5 + nParams];
0955   edm::LogInfo("Vx3DHLTAnalyzer") << "mean x --> " << fitResults[6] << " +/- " << fitResults[6 + nParams];
0956   edm::LogInfo("Vx3DHLTAnalyzer") << "mean y --> " << fitResults[7] << " +/- " << fitResults[7 + nParams];
0957   edm::LogInfo("Vx3DHLTAnalyzer") << "mean z --> " << fitResults[8] << " +/- " << fitResults[8 + nParams];
0958 }
0959 
0960 void Vx3DHLTAnalyzer::dqmBeginLuminosityBlock(const LuminosityBlock& lumiBlock, const EventSetup& iSetup) {
0961   // @@@ If statement to avoid problems with non-sequential lumisections @@@
0962   if ((lumiCounter == 0) && (lumiBlock.luminosityBlock() > lastLumiOfFit)) {
0963     beginTimeOfFit = lumiBlock.beginTime().value();
0964     beginLumiOfFit = lumiBlock.luminosityBlock();
0965     lumiCounter++;
0966   } else if ((lumiCounter != 0) && (lumiBlock.luminosityBlock() >= (beginLumiOfFit + lumiCounter)))
0967     lumiCounter++;
0968   else
0969     reset("scratch");
0970 }
0971 
0972 void Vx3DHLTAnalyzer::dqmEndLuminosityBlock(const LuminosityBlock& lumiBlock, const EventSetup& iSetup) {
0973   stringstream histTitle;
0974   double minXfit, maxXfit;
0975   int goodData;
0976 
0977   if ((nLumiFit != 0) && (lumiCounter % nLumiFit == 0) && (beginTimeOfFit != 0) && (runNumber != 0)) {
0978     endTimeOfFit = lumiBlock.endTime().value();
0979     endLumiOfFit = lumiBlock.luminosityBlock();
0980     lastLumiOfFit = endLumiOfFit;
0981     vector<double> vals;
0982 
0983     hitCounter->getTH1()->SetBinContent(lastLumiOfFit, (double)totalHits);
0984     hitCounter->getTH1()->SetBinError(
0985         lastLumiOfFit,
0986         (totalHits != 0 ? 1.
0987                         : 0.));  // It's not sqrt(n) because we want to weight all entries in the same way for the fit
0988 
0989     if (dataFromFit == true) {
0990       vector<double> fitResults;
0991 
0992       fitResults.push_back(Vx_X->getTH1()->GetRMS() * Vx_X->getTH1()->GetRMS());
0993       fitResults.push_back(Vx_Y->getTH1()->GetRMS() * Vx_Y->getTH1()->GetRMS());
0994       fitResults.push_back(Vx_Z->getTH1()->GetRMS() * Vx_Z->getTH1()->GetRMS());
0995       fitResults.push_back(0.0);
0996       fitResults.push_back(0.0);
0997       fitResults.push_back(0.0);
0998       fitResults.push_back(Vx_X->getTH1()->GetMean());
0999       fitResults.push_back(Vx_Y->getTH1()->GetMean());
1000       fitResults.push_back(Vx_Z->getTH1()->GetMean());
1001       for (unsigned int i = 0; i < nParams; i++)
1002         fitResults.push_back(0.0);
1003 
1004       if (internalDebug == true) {
1005         edm::LogInfo("Vx3DHLTAnalyzer") << "\t@@@ Beam Spot parameters - prefit @@@";
1006 
1007         printFitParams(fitResults);
1008 
1009         edm::LogInfo("Vx3DHLTAnalyzer") << "Runnumber " << runNumber;
1010         edm::LogInfo("Vx3DHLTAnalyzer") << "BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) << " "
1011                                         << (beginTimeOfFit >> 32);
1012         edm::LogInfo("Vx3DHLTAnalyzer") << "EndTimeOfFit " << formatTime(endTimeOfFit >> 32) << " "
1013                                         << (endTimeOfFit >> 32);
1014         edm::LogInfo("Vx3DHLTAnalyzer") << "LumiRange " << beginLumiOfFit << " - " << endLumiOfFit;
1015       }
1016 
1017       goodData = MyFit(&fitResults);
1018 
1019       if (internalDebug == true) {
1020         edm::LogInfo("Vx3DHLTAnalyzer") << "\t@@@ Beam Spot parameters - postfit @@@";
1021 
1022         printFitParams(fitResults);
1023 
1024         edm::LogInfo("Vx3DHLTAnalyzer") << "goodData --> " << goodData;
1025         edm::LogInfo("Vx3DHLTAnalyzer") << "Used vertices --> " << counterVx;
1026       }
1027 
1028       if (goodData == 0) {
1029         vals.push_back(fitResults[6]);
1030         vals.push_back(fitResults[7]);
1031         vals.push_back(fitResults[8]);
1032         vals.push_back(std::sqrt(std::fabs(fitResults[0])));
1033         vals.push_back(std::sqrt(std::fabs(fitResults[1])));
1034         vals.push_back(std::sqrt(std::fabs(fitResults[2])));
1035         vals.push_back(fitResults[5]);
1036         vals.push_back(fitResults[4]);
1037 
1038         vals.push_back(std::pow(fitResults[6 + nParams], 2.));
1039         vals.push_back(std::pow(fitResults[7 + nParams], 2.));
1040         vals.push_back(std::pow(fitResults[8 + nParams], 2.));
1041         vals.push_back(std::pow(std::fabs(fitResults[0 + nParams]) / (2. * std::sqrt(std::fabs(fitResults[0]))), 2.));
1042         vals.push_back(std::pow(std::fabs(fitResults[1 + nParams]) / (2. * std::sqrt(std::fabs(fitResults[1]))), 2.));
1043         vals.push_back(std::pow(std::fabs(fitResults[2 + nParams]) / (2. * std::sqrt(std::fabs(fitResults[2]))), 2.));
1044         vals.push_back(std::pow(fitResults[5 + nParams], 2.));
1045         vals.push_back(std::pow(fitResults[4 + nParams], 2.));
1046       } else
1047         for (unsigned int i = 0; i < (nParams - 1) * 2; i++)
1048           vals.push_back(0.0);
1049 
1050       fitResults.clear();
1051     } else {
1052       counterVx = Vx_X->getTH1F()->GetEntries();
1053 
1054       if (Vx_X->getTH1F()->GetEntries() >= minNentries) {
1055         goodData = 0;
1056 
1057         vals.push_back(Vx_X->getTH1F()->GetMean());
1058         vals.push_back(Vx_Y->getTH1F()->GetMean());
1059         vals.push_back(Vx_Z->getTH1F()->GetMean());
1060         vals.push_back(Vx_X->getTH1F()->GetRMS());
1061         vals.push_back(Vx_Y->getTH1F()->GetRMS());
1062         vals.push_back(Vx_Z->getTH1F()->GetRMS());
1063         vals.push_back(0.0);
1064         vals.push_back(0.0);
1065 
1066         vals.push_back(std::pow(Vx_X->getTH1F()->GetMeanError(), 2.));
1067         vals.push_back(std::pow(Vx_Y->getTH1F()->GetMeanError(), 2.));
1068         vals.push_back(std::pow(Vx_Z->getTH1F()->GetMeanError(), 2.));
1069         vals.push_back(std::pow(Vx_X->getTH1F()->GetRMSError(), 2.));
1070         vals.push_back(std::pow(Vx_Y->getTH1F()->GetRMSError(), 2.));
1071         vals.push_back(std::pow(Vx_Z->getTH1F()->GetRMSError(), 2.));
1072         vals.push_back(0.0);
1073         vals.push_back(0.0);
1074       } else {
1075         goodData = -2;
1076         for (unsigned int i = 0; i < (nParams - 1) * 2; i++)
1077           vals.push_back(0.0);
1078       }
1079     }
1080 
1081     // vals[0]  = X0
1082     // vals[1]  = Y0
1083     // vals[2]  = Z0
1084     // vals[3]  = sigmaX0
1085     // vals[4]  = sigmaY0
1086     // vals[5]  = sigmaZ0
1087     // vals[6]  = dxdz
1088     // vals[7]  = dydz
1089 
1090     // vals[8]  = err^2 X0
1091     // vals[9]  = err^2 Y0
1092     // vals[10] = err^2 Z0
1093     // vals[11] = err^2 sigmaX0
1094     // vals[12] = err^2 sigmaY0
1095     // vals[13] = err^2 sigmaZ0
1096     // vals[14] = err^2 dxdz
1097     // vals[15] = err^2 dydz
1098 
1099     numberFits++;
1100     writeToFile(&vals, beginTimeOfFit, endTimeOfFit, beginLumiOfFit, endLumiOfFit, 3);
1101     if (internalDebug == true)
1102       edm::LogInfo("Vx3DHLTAnalyzer") << "\tUsed vertices: " << counterVx;
1103 
1104     statusCounter->getTH1()->SetBinContent(lastLumiOfFit, (double)goodData);
1105     statusCounter->getTH1()->SetBinError(lastLumiOfFit, 1e-3);
1106 
1107     // Copy vertex position histograms into to-fit histograms
1108     if (goodData == 0)
1109       reset("fit");
1110     else if (lumiCounter >= maxLumiIntegration) {
1111       reset("fit");
1112       reset("whole");
1113     }
1114 
1115     for (int i = 0; i < Vx_X_Fit->getTH1()->GetNbinsX(); i++) {
1116       Vx_X_Fit->getTH1()->SetBinContent(
1117           i + 1, Vx_X_Fit->getTH1()->GetBinContent(i + 1) + Vx_X->getTH1()->GetBinContent(i + 1));
1118       Vx_X_Fit->getTH1()->SetBinError(i + 1, sqrt(Vx_X_Fit->getTH1()->GetBinContent(i + 1)));
1119     }
1120 
1121     for (int i = 0; i < Vx_Y_Fit->getTH1()->GetNbinsX(); i++) {
1122       Vx_Y_Fit->getTH1()->SetBinContent(
1123           i + 1, Vx_Y_Fit->getTH1()->GetBinContent(i + 1) + Vx_Y->getTH1()->GetBinContent(i + 1));
1124       Vx_Y_Fit->getTH1()->SetBinError(i + 1, sqrt(Vx_Y_Fit->getTH1()->GetBinContent(i + 1)));
1125     }
1126 
1127     for (int i = 0; i < Vx_Z_Fit->getTH1()->GetNbinsX(); i++) {
1128       Vx_Z_Fit->getTH1()->SetBinContent(
1129           i + 1, Vx_Z_Fit->getTH1()->GetBinContent(i + 1) + Vx_Z->getTH1()->GetBinContent(i + 1));
1130       Vx_Z_Fit->getTH1()->SetBinError(i + 1, sqrt(Vx_Z_Fit->getTH1()->GetBinContent(i + 1)));
1131     }
1132 
1133     // Check data quality
1134     if (goodData == 0) {
1135       numberGoodFits++;
1136 
1137       histTitle << "Ongoing: fitted lumis " << beginLumiOfFit << " - " << endLumiOfFit;
1138       reset("whole");
1139     } else {
1140       if (goodData == -2)
1141         histTitle << "Ongoing: not enough evts (" << lumiCounter << " - " << maxLumiIntegration << " lumis)";
1142       else
1143         histTitle << "Ongoing: temporary problems (" << lumiCounter << " - " << maxLumiIntegration << " lumis)";
1144 
1145       if (lumiCounter >= maxLumiIntegration) {
1146         statusCounter->getTH1()->SetBinContent(lastLumiOfFit, -5);
1147         statusCounter->getTH1()->SetBinError(lastLumiOfFit, 1e-3);
1148       } else
1149         reset("hitCounter");
1150     }
1151 
1152     reportSummary->Fill((numberFits != 0 ? ((double)numberGoodFits) / ((double)numberFits) : -1));
1153     reportSummaryMap->getTH1()->SetBinContent(
1154         1, 1, (numberFits != 0 ? ((double)numberGoodFits) / ((double)numberFits) : -1));
1155 
1156     fitResults->setAxisTitle(histTitle.str(), 1);
1157 
1158     fitResults->setBinContent(1, 9, vals[0]);
1159     fitResults->setBinContent(1, 8, vals[1]);
1160     fitResults->setBinContent(1, 7, vals[2]);
1161     fitResults->setBinContent(1, 6, vals[3]);
1162     fitResults->setBinContent(1, 5, vals[4]);
1163     fitResults->setBinContent(1, 4, vals[5]);
1164     fitResults->setBinContent(1, 3, vals[6]);
1165     fitResults->setBinContent(1, 2, vals[7]);
1166     fitResults->setBinContent(1, 1, counterVx);
1167 
1168     fitResults->setBinContent(2, 9, std::sqrt(vals[8]));
1169     fitResults->setBinContent(2, 8, std::sqrt(vals[9]));
1170     fitResults->setBinContent(2, 7, std::sqrt(vals[10]));
1171     fitResults->setBinContent(2, 6, std::sqrt(vals[11]));
1172     fitResults->setBinContent(2, 5, std::sqrt(vals[12]));
1173     fitResults->setBinContent(2, 4, std::sqrt(vals[13]));
1174     fitResults->setBinContent(2, 3, std::sqrt(vals[14]));
1175     fitResults->setBinContent(2, 2, std::sqrt(vals[15]));
1176     fitResults->setBinContent(2, 1, std::sqrt(counterVx));
1177 
1178     // Linear fit to the historical plots
1179     TF1* myLinFit = new TF1(
1180         "myLinFit", "[0] + [1]*x", mXlumi->getTH1()->GetXaxis()->GetXmin(), mXlumi->getTH1()->GetXaxis()->GetXmax());
1181     myLinFit->SetLineColor(2);
1182     myLinFit->SetLineWidth(2);
1183     myLinFit->SetParName(0, "Inter.");
1184     myLinFit->SetParName(1, "Slope");
1185 
1186     mXlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[0]);
1187     mXlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[8]));
1188     myLinFit->SetParameter(0, mXlumi->getTH1()->GetMean(2));
1189     myLinFit->SetParameter(1, 0.0);
1190     mXlumi->getTH1()->Fit(myLinFit, "QR");
1191 
1192     mYlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[1]);
1193     mYlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[9]));
1194     myLinFit->SetParameter(0, mYlumi->getTH1()->GetMean(2));
1195     myLinFit->SetParameter(1, 0.0);
1196     mYlumi->getTH1()->Fit(myLinFit, "QR");
1197 
1198     mZlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[2]);
1199     mZlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[10]));
1200     myLinFit->SetParameter(0, mZlumi->getTH1()->GetMean(2));
1201     myLinFit->SetParameter(1, 0.0);
1202     mZlumi->getTH1()->Fit(myLinFit, "QR");
1203 
1204     sXlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[3]);
1205     sXlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[11]));
1206     myLinFit->SetParameter(0, sXlumi->getTH1()->GetMean(2));
1207     myLinFit->SetParameter(1, 0.0);
1208     sXlumi->getTH1()->Fit(myLinFit, "QR");
1209 
1210     sYlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[4]);
1211     sYlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[12]));
1212     myLinFit->SetParameter(0, sYlumi->getTH1()->GetMean(2));
1213     myLinFit->SetParameter(1, 0.0);
1214     sYlumi->getTH1()->Fit(myLinFit, "QR");
1215 
1216     sZlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[5]);
1217     sZlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[13]));
1218     myLinFit->SetParameter(0, sZlumi->getTH1()->GetMean(2));
1219     myLinFit->SetParameter(1, 0.0);
1220     sZlumi->getTH1()->Fit(myLinFit, "QR");
1221 
1222     dxdzlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[6]);
1223     dxdzlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[14]));
1224     myLinFit->SetParameter(0, dxdzlumi->getTH1()->GetMean(2));
1225     myLinFit->SetParameter(1, 0.0);
1226     dxdzlumi->getTH1()->Fit(myLinFit, "QR");
1227 
1228     dydzlumi->getTH1()->SetBinContent(lastLumiOfFit, vals[7]);
1229     dydzlumi->getTH1()->SetBinError(lastLumiOfFit, std::sqrt(vals[15]));
1230     myLinFit->SetParameter(0, dydzlumi->getTH1()->GetMean(2));
1231     myLinFit->SetParameter(1, 0.0);
1232     dydzlumi->getTH1()->Fit(myLinFit, "QR");
1233 
1234     myLinFit->SetParameter(0, hitCounter->getTH1()->GetMean(2));
1235     myLinFit->SetParameter(1, 0.0);
1236     hitCounter->getTH1()->Fit(myLinFit, "QR");
1237 
1238     goodVxCounter->getTH1()->SetBinContent(lastLumiOfFit, (double)counterVx);
1239     goodVxCounter->getTH1()->SetBinError(
1240         lastLumiOfFit,
1241         (counterVx != 0 ? 1.
1242                         : 0.));  // It's not sqrt(n) because we want to weight all entries in the same way for the fit
1243     myLinFit->SetParameter(0, goodVxCounter->getTH1()->GetMean(2));
1244     myLinFit->SetParameter(1, 0.0);
1245     goodVxCounter->getTH1()->Fit(myLinFit, "QR");
1246 
1247     delete myLinFit;
1248     vals.clear();
1249 
1250     // Gaussian fit to 1D vertex coordinate distributions
1251     TF1* myGaussFit = new TF1("myGaussFit",
1252                               "[0]*exp(-(x-[1])*(x-[1])/(2*[2]*[2]))",
1253                               Vx_Z_Fit->getTH1()->GetXaxis()->GetXmin(),
1254                               Vx_Z_Fit->getTH1()->GetXaxis()->GetXmax());
1255     myGaussFit->SetLineColor(2);
1256     myGaussFit->SetLineWidth(2);
1257     myGaussFit->SetParName(0, "Ampl.");
1258     myGaussFit->SetParName(1, "#mu");
1259     myGaussFit->SetParName(2, "#sigma");
1260 
1261     myGaussFit->SetParameter(0, Vx_X_Fit->getTH1()->GetMaximum());
1262     myGaussFit->SetParameter(1, Vx_X_Fit->getTH1()->GetMean());
1263     myGaussFit->SetParameter(2, Vx_X_Fit->getTH1()->GetRMS());
1264     minXfit = Vx_X_Fit->getTH1()->GetBinLowEdge(1);
1265     for (int i = 0; i < Vx_X_Fit->getTH1()->GetNbinsX(); i++) {
1266       if (Vx_X_Fit->getTH1()->GetBinContent(i + 1) > 0) {
1267         minXfit = Vx_X_Fit->getTH1()->GetBinLowEdge(i + 1);
1268         break;
1269       }
1270     }
1271     maxXfit = Vx_X_Fit->getTH1()->GetBinLowEdge(Vx_X_Fit->getTH1()->GetNbinsX());
1272     for (int i = Vx_X_Fit->getTH1()->GetNbinsX(); i > 0; i--) {
1273       if (Vx_X_Fit->getTH1()->GetBinContent(i) > 0) {
1274         maxXfit = Vx_X_Fit->getTH1()->GetBinLowEdge(i);
1275         break;
1276       }
1277     }
1278     myGaussFit->SetRange(minXfit - (maxXfit - minXfit) / 2., maxXfit + (maxXfit - minXfit) / 2.);
1279     if (Vx_X_Fit->getTH1()->GetEntries() > 0)
1280       Vx_X_Fit->getTH1()->Fit(myGaussFit, "QRL");
1281 
1282     myGaussFit->SetParameter(0, Vx_Y_Fit->getTH1()->GetMaximum());
1283     myGaussFit->SetParameter(1, Vx_Y_Fit->getTH1()->GetMean());
1284     myGaussFit->SetParameter(2, Vx_Y_Fit->getTH1()->GetRMS());
1285     minXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(1);
1286     for (int i = 0; i < Vx_Y_Fit->getTH1()->GetNbinsX(); i++) {
1287       if (Vx_Y_Fit->getTH1()->GetBinContent(i + 1) > 0) {
1288         minXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(i + 1);
1289         break;
1290       }
1291     }
1292     maxXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(Vx_Y_Fit->getTH1()->GetNbinsX());
1293     for (int i = Vx_Y_Fit->getTH1()->GetNbinsX(); i > 0; i--) {
1294       if (Vx_Y_Fit->getTH1()->GetBinContent(i) > 0) {
1295         maxXfit = Vx_Y_Fit->getTH1()->GetBinLowEdge(i);
1296         break;
1297       }
1298     }
1299     myGaussFit->SetRange(minXfit - (maxXfit - minXfit) / 2., maxXfit + (maxXfit - minXfit) / 2.);
1300     if (Vx_Y_Fit->getTH1()->GetEntries() > 0)
1301       Vx_Y_Fit->getTH1()->Fit(myGaussFit, "QRL");
1302 
1303     myGaussFit->SetParameter(0, Vx_Z_Fit->getTH1()->GetMaximum());
1304     myGaussFit->SetParameter(1, Vx_Z_Fit->getTH1()->GetMean());
1305     myGaussFit->SetParameter(2, Vx_Z_Fit->getTH1()->GetRMS());
1306     minXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(1);
1307     for (int i = 0; i < Vx_Z_Fit->getTH1()->GetNbinsX(); i++) {
1308       if (Vx_Z_Fit->getTH1()->GetBinContent(i + 1) > 0) {
1309         minXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(i + 1);
1310         break;
1311       }
1312     }
1313     maxXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(Vx_Z_Fit->getTH1()->GetNbinsX());
1314     for (int i = Vx_Z_Fit->getTH1()->GetNbinsX(); i > 0; i--) {
1315       if (Vx_Z_Fit->getTH1()->GetBinContent(i) > 0) {
1316         maxXfit = Vx_Z_Fit->getTH1()->GetBinLowEdge(i);
1317         break;
1318       }
1319     }
1320     myGaussFit->SetRange(minXfit - (maxXfit - minXfit) / 2., maxXfit + (maxXfit - minXfit) / 2.);
1321     if (Vx_Z_Fit->getTH1()->GetEntries() > 0)
1322       Vx_Z_Fit->getTH1()->Fit(myGaussFit, "QRL");
1323 
1324     delete myGaussFit;
1325   } else if ((nLumiFit != 0) && (lumiCounter % nLumiFit != 0) && (beginTimeOfFit != 0) && (runNumber != 0)) {
1326     histTitle << "Ongoing: accumulating evts (" << lumiCounter % nLumiFit << " - " << nLumiFit << " in " << lumiCounter
1327               << " - " << maxLumiIntegration << " lumis)";
1328     fitResults->setAxisTitle(histTitle.str(), 1);
1329     if ((debugMode == true) && (outputDebugFile.is_open() == true)) {
1330       outputDebugFile << "\n"
1331                       << "Runnumber " << runNumber << endl;
1332       outputDebugFile << "BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) << " " << (beginTimeOfFit >> 32) << endl;
1333       outputDebugFile << "BeginLumiRange " << beginLumiOfFit << endl;
1334       outputDebugFile << histTitle.str().c_str() << "\n" << endl;
1335     }
1336   } else if ((nLumiFit == 0) || (beginTimeOfFit == 0) || (runNumber == 0)) {
1337     histTitle << "Ongoing: no ongoing fits";
1338     fitResults->setAxisTitle(histTitle.str(), 1);
1339     if ((debugMode == true) && (outputDebugFile.is_open() == true))
1340       outputDebugFile << histTitle.str().c_str() << "\n" << endl;
1341 
1342     endLumiOfFit = lumiBlock.luminosityBlock();
1343 
1344     hitCounter->getTH1()->SetBinContent(endLumiOfFit, (double)totalHits);
1345     hitCounter->getTH1()->SetBinError(endLumiOfFit, std::sqrt((double)totalHits));
1346 
1347     reset("whole");
1348   }
1349 
1350   if (internalDebug == true)
1351     edm::LogInfo("Vx3DHLTAnalyzer") << "::\tHistogram title: " << histTitle.str();
1352 }
1353 
1354 void Vx3DHLTAnalyzer::bookHistograms(DQMStore::IBooker& ibooker, Run const& iRun, EventSetup const& /* iSetup */) {
1355   ibooker.setCurrentFolder("BeamPixel");
1356 
1357   Vx_X = ibooker.book1D(
1358       "F - vertex x", "Primary Vertex X Distribution", int(rint(xRange / xStep)), -xRange / 2., xRange / 2.);
1359   Vx_Y = ibooker.book1D(
1360       "F - vertex y", "Primary Vertex Y Distribution", int(rint(yRange / yStep)), -yRange / 2., yRange / 2.);
1361   Vx_Z = ibooker.book1D(
1362       "F - vertex z", "Primary Vertex Z Distribution", int(rint(zRange / zStep)), -zRange / 2., zRange / 2.);
1363   Vx_X->setAxisTitle("Primary Vertices X [cm]", 1);
1364   Vx_X->setAxisTitle("Entries [#]", 2);
1365   Vx_Y->setAxisTitle("Primary Vertices Y [cm]", 1);
1366   Vx_Y->setAxisTitle("Entries [#]", 2);
1367   Vx_Z->setAxisTitle("Primary Vertices Z [cm]", 1);
1368   Vx_Z->setAxisTitle("Entries [#]", 2);
1369 
1370   Vx_X_Fit = ibooker.book1D("G - vertex x fit",
1371                             "Primary Vertex X Distribution (For Fit)",
1372                             int(rint(xRange / xStep)),
1373                             -xRange / 2.,
1374                             xRange / 2.);
1375   Vx_Y_Fit = ibooker.book1D("G - vertex y fit",
1376                             "Primary Vertex Y Distribution (For Fit)",
1377                             int(rint(yRange / yStep)),
1378                             -yRange / 2.,
1379                             yRange / 2.);
1380   Vx_Z_Fit = ibooker.book1D("G - vertex z fit",
1381                             "Primary Vertex Z Distribution (For Fit)",
1382                             int(rint(zRange / zStep)),
1383                             -zRange / 2.,
1384                             zRange / 2.);
1385   Vx_X_Fit->setAxisTitle("Primary Vertices X [cm]", 1);
1386   Vx_X_Fit->setAxisTitle("Entries [#]", 2);
1387   Vx_Y_Fit->setAxisTitle("Primary Vertices Y [cm]", 1);
1388   Vx_Y_Fit->setAxisTitle("Entries [#]", 2);
1389   Vx_Z_Fit->setAxisTitle("Primary Vertices Z [cm]", 1);
1390   Vx_Z_Fit->setAxisTitle("Entries [#]", 2);
1391 
1392   Vx_X_Cum = ibooker.book1D("I - vertex x cum",
1393                             "Primary Vertex X Distribution (Cumulative)",
1394                             int(rint(xRange / xStep)),
1395                             -xRange / 2.,
1396                             xRange / 2.);
1397   Vx_Y_Cum = ibooker.book1D("I - vertex y cum",
1398                             "Primary Vertex Y Distribution (Cumulative)",
1399                             int(rint(yRange / yStep)),
1400                             -yRange / 2.,
1401                             yRange / 2.);
1402   Vx_Z_Cum = ibooker.book1D("I - vertex z cum",
1403                             "Primary Vertex Z Distribution (Cumulative)",
1404                             int(rint(zRange / zStep)),
1405                             -zRange / 2.,
1406                             zRange / 2.);
1407   Vx_X_Cum->setAxisTitle("Primary Vertices X [cm]", 1);
1408   Vx_X_Cum->setAxisTitle("Entries [#]", 2);
1409   Vx_Y_Cum->setAxisTitle("Primary Vertices Y [cm]", 1);
1410   Vx_Y_Cum->setAxisTitle("Entries [#]", 2);
1411   Vx_Z_Cum->setAxisTitle("Primary Vertices Z [cm]", 1);
1412   Vx_Z_Cum->setAxisTitle("Entries [#]", 2);
1413 
1414   mXlumi = ibooker.book1D(
1415       "B - muX vs lumi", "#mu_{x} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1416   mYlumi = ibooker.book1D(
1417       "B - muY vs lumi", "#mu_{y} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1418   mZlumi = ibooker.book1D(
1419       "B - muZ vs lumi", "#mu_{z} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1420   mXlumi->setAxisTitle("Lumisection [#]", 1);
1421   mXlumi->setAxisTitle("#mu_{x} [cm]", 2);
1422   mXlumi->getTH1()->SetOption("E1");
1423   mYlumi->setAxisTitle("Lumisection [#]", 1);
1424   mYlumi->setAxisTitle("#mu_{y} [cm]", 2);
1425   mYlumi->getTH1()->SetOption("E1");
1426   mZlumi->setAxisTitle("Lumisection [#]", 1);
1427   mZlumi->setAxisTitle("#mu_{z} [cm]", 2);
1428   mZlumi->getTH1()->SetOption("E1");
1429 
1430   sXlumi = ibooker.book1D(
1431       "C - sigmaX vs lumi", "#sigma_{x} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1432   sYlumi = ibooker.book1D(
1433       "C - sigmaY vs lumi", "#sigma_{y} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1434   sZlumi = ibooker.book1D(
1435       "C - sigmaZ vs lumi", "#sigma_{z} vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1436   sXlumi->setAxisTitle("Lumisection [#]", 1);
1437   sXlumi->setAxisTitle("#sigma_{x} [cm]", 2);
1438   sXlumi->getTH1()->SetOption("E1");
1439   sYlumi->setAxisTitle("Lumisection [#]", 1);
1440   sYlumi->setAxisTitle("#sigma_{y} [cm]", 2);
1441   sYlumi->getTH1()->SetOption("E1");
1442   sZlumi->setAxisTitle("Lumisection [#]", 1);
1443   sZlumi->setAxisTitle("#sigma_{z} [cm]", 2);
1444   sZlumi->getTH1()->SetOption("E1");
1445 
1446   dxdzlumi = ibooker.book1D(
1447       "D - dxdz vs lumi", "dX/dZ vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1448   dydzlumi = ibooker.book1D(
1449       "D - dydz vs lumi", "dY/dZ vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1450   dxdzlumi->setAxisTitle("Lumisection [#]", 1);
1451   dxdzlumi->setAxisTitle("dX/dZ [rad]", 2);
1452   dxdzlumi->getTH1()->SetOption("E1");
1453   dydzlumi->setAxisTitle("Lumisection [#]", 1);
1454   dydzlumi->setAxisTitle("dY/dZ [rad]", 2);
1455   dydzlumi->getTH1()->SetOption("E1");
1456 
1457   Vx_ZX = ibooker.book2D("E - vertex zx",
1458                          "Primary Vertex ZX Distribution",
1459                          int(rint(zRange / zStep)),
1460                          -zRange / 2.,
1461                          zRange / 2.,
1462                          int(rint(xRange / xStep)),
1463                          -xRange / 2.,
1464                          xRange / 2.);
1465   Vx_ZY = ibooker.book2D("E - vertex zy",
1466                          "Primary Vertex ZY Distribution",
1467                          int(rint(zRange / zStep)),
1468                          -zRange / 2.,
1469                          zRange / 2.,
1470                          int(rint(yRange / yStep)),
1471                          -yRange / 2.,
1472                          yRange / 2.);
1473   Vx_XY = ibooker.book2D("E - vertex xy",
1474                          "Primary Vertex XY Distribution",
1475                          int(rint(xRange / xStep)),
1476                          -xRange / 2.,
1477                          xRange / 2.,
1478                          int(rint(yRange / yStep)),
1479                          -yRange / 2.,
1480                          yRange / 2.);
1481   Vx_ZX->setAxisTitle("Primary Vertices Z [cm]", 1);
1482   Vx_ZX->setAxisTitle("Primary Vertices X [cm]", 2);
1483   Vx_ZX->setAxisTitle("Entries [#]", 3);
1484   Vx_ZY->setAxisTitle("Primary Vertices Z [cm]", 1);
1485   Vx_ZY->setAxisTitle("Primary Vertices Y [cm]", 2);
1486   Vx_ZY->setAxisTitle("Entries [#]", 3);
1487   Vx_XY->setAxisTitle("Primary Vertices X [cm]", 1);
1488   Vx_XY->setAxisTitle("Primary Vertices Y [cm]", 2);
1489   Vx_XY->setAxisTitle("Entries [#]", 3);
1490 
1491   Vx_ZX_Cum = ibooker.book2D("H - vertex zx cum",
1492                              "Primary Vertex ZX Distribution (Cumulative)",
1493                              int(rint(zRange / zStep)),
1494                              -zRange / 2.,
1495                              zRange / 2.,
1496                              int(rint(xRange / xStep)),
1497                              -xRange / 2.,
1498                              xRange / 2.);
1499   Vx_ZY_Cum = ibooker.book2D("H - vertex zy cum",
1500                              "Primary Vertex ZY Distribution (Cumulative)",
1501                              int(rint(zRange / zStep)),
1502                              -zRange / 2.,
1503                              zRange / 2.,
1504                              int(rint(yRange / yStep)),
1505                              -yRange / 2.,
1506                              yRange / 2.);
1507   Vx_XY_Cum = ibooker.book2D("H - vertex xy cum",
1508                              "Primary Vertex XY Distribution (Cumulative)",
1509                              int(rint(xRange / xStep)),
1510                              -xRange / 2.,
1511                              xRange / 2.,
1512                              int(rint(yRange / yStep)),
1513                              -yRange / 2.,
1514                              yRange / 2.);
1515   Vx_ZX_Cum->setAxisTitle("Primary Vertices Z [cm]", 1);
1516   Vx_ZX_Cum->setAxisTitle("Primary Vertices X [cm]", 2);
1517   Vx_ZX_Cum->setAxisTitle("Entries [#]", 3);
1518   Vx_ZY_Cum->setAxisTitle("Primary Vertices Z [cm]", 1);
1519   Vx_ZY_Cum->setAxisTitle("Primary Vertices Y [cm]", 2);
1520   Vx_ZY_Cum->setAxisTitle("Entries [#]", 3);
1521   Vx_XY_Cum->setAxisTitle("Primary Vertices X [cm]", 1);
1522   Vx_XY_Cum->setAxisTitle("Primary Vertices Y [cm]", 2);
1523   Vx_XY_Cum->setAxisTitle("Entries [#]", 3);
1524 
1525   hitCounter = ibooker.book1D(
1526       "J - pixelHits vs lumi", "# Pixel-Hits vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1527   hitCounter->setAxisTitle("Lumisection [#]", 1);
1528   hitCounter->setAxisTitle("Pixel-Hits [#]", 2);
1529   hitCounter->getTH1()->SetOption("E1");
1530 
1531   goodVxCounter = ibooker.book1D("K - good vertices vs lumi",
1532                                  "# Good vertices vs. Lumisection",
1533                                  nLumiXaxisRange,
1534                                  0.5,
1535                                  ((double)nLumiXaxisRange) + 0.5);
1536   goodVxCounter->setAxisTitle("Lumisection [#]", 1);
1537   goodVxCounter->setAxisTitle("Good vertices [#]", 2);
1538   goodVxCounter->getTH1()->SetOption("E1");
1539 
1540   statusCounter = ibooker.book1D(
1541       "L - status vs lumi", "App. Status vs. Lumisection", nLumiXaxisRange, 0.5, ((double)nLumiXaxisRange) + 0.5);
1542   statusCounter->setAxisTitle("Lumisection [#]", 1);
1543   statusCounter->getTH1()->SetOption("E1");
1544   statusCounter->getTH1()->GetYaxis()->Set(12, -6.5, 5.5);
1545   statusCounter->getTH1()->GetYaxis()->SetBinLabel(1, "FatalRootError");
1546   statusCounter->getTH1()->GetYaxis()->SetBinLabel(2, "Max Lumi.");
1547   statusCounter->getTH1()->GetYaxis()->SetBinLabel(3, "Neg. det.");
1548   statusCounter->getTH1()->GetYaxis()->SetBinLabel(4, "Infinite err.");
1549   statusCounter->getTH1()->GetYaxis()->SetBinLabel(5, "No vtx.");
1550   statusCounter->getTH1()->GetYaxis()->SetBinLabel(6, "Infinite EDM");
1551   statusCounter->getTH1()->GetYaxis()->SetBinLabel(7, "OK");
1552   statusCounter->getTH1()->GetYaxis()->SetBinLabel(8, "MINUIT stat.");
1553   statusCounter->getTH1()->GetYaxis()->SetBinLabel(9, "MINUIT stat.");
1554   statusCounter->getTH1()->GetYaxis()->SetBinLabel(10, "MINUIT stat.");
1555   statusCounter->getTH1()->GetYaxis()->SetBinLabel(11, "MINUIT stat.");
1556   statusCounter->getTH1()->GetYaxis()->SetBinLabel(12, "MINUIT stat.");
1557 
1558   fitResults = ibooker.book2D("A - fit results", "Results of Beam Spot Fit", 2, 0., 2., 9, 0., 9.);
1559   fitResults->setAxisTitle("Ongoing: bootstrapping", 1);
1560   fitResults->setBinLabel(9, "X[cm]", 2);
1561   fitResults->setBinLabel(8, "Y[cm]", 2);
1562   fitResults->setBinLabel(7, "Z[cm]", 2);
1563   fitResults->setBinLabel(6, "#sigma_{X}[cm]", 2);
1564   fitResults->setBinLabel(5, "#sigma_{Y}[cm]", 2);
1565   fitResults->setBinLabel(4, "#sigma_{Z}[cm]", 2);
1566   fitResults->setBinLabel(3, "#frac{dX}{dZ}[rad]", 2);
1567   fitResults->setBinLabel(2, "#frac{dY}{dZ}[rad]", 2);
1568   fitResults->setBinLabel(1, "Vtx[#]", 2);
1569   fitResults->setBinLabel(1, "Value", 1);
1570   fitResults->setBinLabel(2, "Error (stat)", 1);
1571   fitResults->getTH1()->SetOption("text");
1572 
1573   ibooker.setCurrentFolder("BeamPixel/EventInfo");
1574 
1575   reportSummary = ibooker.bookFloat("reportSummary");
1576   reportSummary->Fill(-1);
1577   reportSummaryMap = ibooker.book2D("reportSummaryMap", "Pixel-Vertices Beam Spot: % Good Fits", 1, 0., 1., 1, 0., 1.);
1578   reportSummaryMap->getTH1()->SetBinContent(1, 1, -1);
1579 
1580   ibooker.setCurrentFolder("BeamPixel/EventInfo/reportSummaryContents");
1581 
1582   // Convention for reportSummary and reportSummaryMap:
1583   // - -1%  at the moment of creation of the histogram (i.e. white histogram)
1584   // - n%  numberGoodFits / numberFits
1585 
1586   reset("scratch");  // Initialize histograms after creation
1587 }
1588 
1589 // Define this as a plug-in
1590 DEFINE_FWK_MODULE(Vx3DHLTAnalyzer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!