Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    Vx3DHLTAnalyzer
00004 // Class:      Vx3DHLTAnalyzer
00005 // 
00013 //
00014 // Original Author:  Mauro Dinardo,28 S-020,+41227673777,
00015 //         Created:  Tue Feb 23 13:15:31 CET 2010
00016 // $Id: Vx3DHLTAnalyzer.cc,v 1.102 2012/03/12 11:52:31 dinardo Exp $
00017 
00018 
00019 #include "DQM/BeamMonitor/plugins/Vx3DHLTAnalyzer.h"
00020 
00021 #include "FWCore/ServiceRegistry/interface/Service.h"
00022 #include "FWCore/Framework/interface/LuminosityBlock.h"
00023 
00024 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
00025 #include "DataFormats/TrackReco/interface/Track.h"
00026 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00027 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00028 #include "DataFormats/VertexReco/interface/Vertex.h"
00029 
00030 #include <TFitterMinuit.h>
00031 
00032 
00033 using namespace std;
00034 using namespace reco;
00035 using namespace edm;
00036 
00037 
00038 Vx3DHLTAnalyzer::Vx3DHLTAnalyzer(const ParameterSet& iConfig)
00039 {
00040   vertexCollection = edm::InputTag("pixelVertices");
00041   debugMode        = true;
00042   nLumiReset       = 1;
00043   dataFromFit      = true;
00044   minNentries      = 35;
00045   xRange           = 2.;
00046   xStep            = 0.001;
00047   yRange           = 2.;
00048   yStep            = 0.001;
00049   zRange           = 30.;
00050   zStep            = 0.05;
00051   VxErrCorr        = 1.5;
00052   fileName         = "BeamPixelResults.txt";
00053 
00054   vertexCollection = iConfig.getParameter<InputTag>("vertexCollection");
00055   debugMode        = iConfig.getParameter<bool>("debugMode");
00056   nLumiReset       = iConfig.getParameter<unsigned int>("nLumiReset");
00057   dataFromFit      = iConfig.getParameter<bool>("dataFromFit");
00058   minNentries      = iConfig.getParameter<unsigned int>("minNentries");
00059   xRange           = iConfig.getParameter<double>("xRange");
00060   xStep            = iConfig.getParameter<double>("xStep");
00061   yRange           = iConfig.getParameter<double>("yRange");
00062   yStep            = iConfig.getParameter<double>("yStep");
00063   zRange           = iConfig.getParameter<double>("zRange");
00064   zStep            = iConfig.getParameter<double>("zStep");
00065   VxErrCorr        = iConfig.getParameter<double>("VxErrCorr");
00066   fileName         = iConfig.getParameter<string>("fileName");
00067 }
00068 
00069 
00070 Vx3DHLTAnalyzer::~Vx3DHLTAnalyzer()
00071 {
00072 }
00073 
00074 
00075 void Vx3DHLTAnalyzer::analyze(const Event& iEvent, const EventSetup& iSetup)
00076 {
00077   Handle<VertexCollection> Vx3DCollection;
00078   iEvent.getByLabel(vertexCollection,Vx3DCollection);
00079 
00080   unsigned int i,j;
00081   double det;
00082   VertexType MyVertex;
00083 
00084   if (runNumber != iEvent.id().run())
00085     {
00086       reset("scratch");
00087       runNumber = iEvent.id().run();
00088 
00089       if (debugMode == true)
00090         {
00091           stringstream debugFile;
00092           string tmp(fileName);
00093 
00094           if (outputDebugFile.is_open() == true) outputDebugFile.close();
00095           tmp.erase(strlen(fileName.c_str())-4,4);
00096           debugFile << tmp.c_str() << "_Run" << iEvent.id().run() << ".txt";
00097           outputDebugFile.open(debugFile.str().c_str(), ios::out);
00098           outputDebugFile.close();
00099           outputDebugFile.open(debugFile.str().c_str(), ios::app);
00100         }
00101 
00102       beginLuminosityBlock(iEvent.getLuminosityBlock(),iSetup);
00103     }
00104   else if (beginTimeOfFit != 0)
00105     {
00106       totalHits += HitCounter(iEvent);
00107 
00108       for (vector<Vertex>::const_iterator it3DVx = Vx3DCollection->begin(); it3DVx != Vx3DCollection->end(); it3DVx++)
00109         {
00110           if ((it3DVx->isValid() == true) &&
00111               (it3DVx->isFake() == false) &&
00112               (it3DVx->ndof() >= minVxDoF))
00113             {
00114               for (i = 0; i < DIM; i++)
00115                 {
00116                   for (j = 0; j < DIM; j++)
00117                     {
00118                       MyVertex.Covariance[i][j] = it3DVx->covariance(i,j);
00119                       if (std::isnan(MyVertex.Covariance[i][j]) == true) break;
00120                     }
00121                   if (j != DIM) break;
00122                 }
00123               det = std::fabs(MyVertex.Covariance[0][0])*(std::fabs(MyVertex.Covariance[1][1])*std::fabs(MyVertex.Covariance[2][2]) - MyVertex.Covariance[1][2]*MyVertex.Covariance[1][2]) -
00124                 MyVertex.Covariance[0][1]*(MyVertex.Covariance[0][1]*std::fabs(MyVertex.Covariance[2][2]) - MyVertex.Covariance[0][2]*MyVertex.Covariance[1][2]) +
00125                 MyVertex.Covariance[0][2]*(MyVertex.Covariance[0][1]*MyVertex.Covariance[1][2] - MyVertex.Covariance[0][2]*std::fabs(MyVertex.Covariance[1][1]));
00126               if ((i == DIM) && (det > 0.))
00127                 {
00128                   MyVertex.x = it3DVx->x();
00129                   MyVertex.y = it3DVx->y();
00130                   MyVertex.z = it3DVx->z();
00131                   Vertices.push_back(MyVertex);
00132                 }
00133               else if (internalDebug == true)
00134                 {
00135                   cout << "Vertex discarded !" << endl;
00136                   for (i = 0; i < DIM; i++)
00137                     for (j = 0; j < DIM; j++)
00138                       cout << "(i,j) --> " << i << "," << j << " --> " << MyVertex.Covariance[i][j] << endl;
00139                 }
00140               
00141               Vx_X->Fill(it3DVx->x());
00142               Vx_Y->Fill(it3DVx->y());
00143               Vx_Z->Fill(it3DVx->z());
00144               
00145               Vx_ZX->Fill(it3DVx->z(), it3DVx->x());
00146               Vx_ZY->Fill(it3DVx->z(), it3DVx->y());
00147               Vx_XY->Fill(it3DVx->x(), it3DVx->y());
00148             }
00149         }
00150     }
00151 }
00152 
00153 
00154 unsigned int Vx3DHLTAnalyzer::HitCounter(const Event& iEvent)
00155 {
00156   edm::Handle<SiPixelRecHitCollection> rechitspixel;
00157   iEvent.getByLabel("siPixelRecHits",rechitspixel);
00158 
00159   unsigned int counter = 0;
00160   
00161   for (SiPixelRecHitCollection::const_iterator j = rechitspixel->begin(); j != rechitspixel->end(); j++)
00162     for (edmNew::DetSet<SiPixelRecHit>::const_iterator h = j->begin(); h != j->end(); h++) counter += h->cluster()->size();       
00163   
00164   return counter;
00165 }
00166 
00167 
00168 char* Vx3DHLTAnalyzer::formatTime (const time_t& t)
00169 {
00170   static char ts[25];
00171   strftime(ts, sizeof(ts), "%Y.%m.%d %H:%M:%S %Z", gmtime(&t));
00172 
00173   return ts;
00174 }
00175 
00176 
00177 void Gauss3DFunc(int& /*npar*/, double* /*gin*/, double& fval, double* par, int /*iflag*/)
00178 {
00179   double K[DIM][DIM]; // Covariance Matrix
00180   double M[DIM][DIM]; // K^-1
00181   double det;
00182   double sumlog = 0.;
00183 
00184 //   par[0] = K(0,0) --> Var[X]
00185 //   par[1] = K(1,1) --> Var[Y]
00186 //   par[2] = K(2,2) --> Var[Z]
00187 //   par[3] = K(0,1) = K(1,0) --> Cov[X,Y]
00188 //   par[4] = K(1,2) = K(2,1) --> Cov[Y,Z] --> dy/dz
00189 //   par[5] = K(0,2) = K(2,0) --> Cov[X,Z] --> dx/dz
00190 //   par[6] = mean x
00191 //   par[7] = mean y
00192 //   par[8] = mean z
00193 
00194   counterVx = 0;
00195   for (unsigned int i = 0; i < Vertices.size(); i++)
00196     {
00197       if ((std::sqrt((Vertices[i].x-xPos)*(Vertices[i].x-xPos) + (Vertices[i].y-yPos)*(Vertices[i].y-yPos)) <= maxTransRadius) &&
00198           (std::fabs(Vertices[i].z-zPos) <= maxLongLength))
00199         {
00200           if (considerVxCovariance == true)
00201             {
00202               K[0][0] = std::fabs(par[0]) + VxErrCorr*VxErrCorr * std::fabs(Vertices[i].Covariance[0][0]);
00203               K[1][1] = std::fabs(par[1]) + VxErrCorr*VxErrCorr * std::fabs(Vertices[i].Covariance[1][1]);
00204               K[2][2] = std::fabs(par[2]) + VxErrCorr*VxErrCorr * std::fabs(Vertices[i].Covariance[2][2]);
00205               K[0][1] = K[1][0] = par[3] + VxErrCorr*VxErrCorr * Vertices[i].Covariance[0][1];
00206               K[1][2] = K[2][1] = par[4]*(std::fabs(par[2])-std::fabs(par[1])) - par[5]*par[3] + VxErrCorr*VxErrCorr * Vertices[i].Covariance[1][2];
00207               K[0][2] = K[2][0] = par[5]*(std::fabs(par[2])-std::fabs(par[0])) - par[4]*par[3] + VxErrCorr*VxErrCorr * Vertices[i].Covariance[0][2];
00208             }
00209           else
00210             {
00211               K[0][0] = std::fabs(par[0]);
00212               K[1][1] = std::fabs(par[1]);
00213               K[2][2] = std::fabs(par[2]);
00214               K[0][1] = K[1][0] = par[3];
00215               K[1][2] = K[2][1] = par[4]*(std::fabs(par[2])-std::fabs(par[1])) - par[5]*par[3];
00216               K[0][2] = K[2][0] = par[5]*(std::fabs(par[2])-std::fabs(par[0])) - par[4]*par[3];
00217             }
00218 
00219           det = K[0][0]*(K[1][1]*K[2][2] - K[1][2]*K[1][2]) -
00220                 K[0][1]*(K[0][1]*K[2][2] - K[0][2]*K[1][2]) +
00221                 K[0][2]*(K[0][1]*K[1][2] - K[0][2]*K[1][1]);
00222 
00223           M[0][0] = (K[1][1]*K[2][2] - K[1][2]*K[1][2]) / det;
00224           M[1][1] = (K[0][0]*K[2][2] - K[0][2]*K[0][2]) / det;
00225           M[2][2] = (K[0][0]*K[1][1] - K[0][1]*K[0][1]) / det;
00226           M[0][1] = M[1][0] = (K[0][2]*K[1][2] - K[0][1]*K[2][2]) / det;
00227           M[1][2] = M[2][1] = (K[0][2]*K[0][1] - K[1][2]*K[0][0]) / det;
00228           M[0][2] = M[2][0] = (K[0][1]*K[1][2] - K[0][2]*K[1][1]) / det;
00229           
00230           sumlog += double(DIM)*std::log(2.*pi) + std::log(std::fabs(det)) +
00231             (M[0][0]*(Vertices[i].x-par[6])*(Vertices[i].x-par[6]) +
00232              M[1][1]*(Vertices[i].y-par[7])*(Vertices[i].y-par[7]) +
00233              M[2][2]*(Vertices[i].z-par[8])*(Vertices[i].z-par[8]) +
00234              2.*M[0][1]*(Vertices[i].x-par[6])*(Vertices[i].y-par[7]) +
00235              2.*M[1][2]*(Vertices[i].y-par[7])*(Vertices[i].z-par[8]) +
00236              2.*M[0][2]*(Vertices[i].x-par[6])*(Vertices[i].z-par[8]));
00237           
00238           counterVx++;
00239         }
00240     }
00241   
00242   fval = sumlog;
00243 }
00244 
00245 
00246 int Vx3DHLTAnalyzer::MyFit(vector<double>* vals)
00247 {
00248   // RETURN CODE:
00249   //  0 == OK
00250   // -2 == NO OK - not enough "minNentries"
00251   // Any other number == NO OK
00252   unsigned int nParams = 9;
00253  
00254   if ((vals != NULL) && (vals->size() == nParams*2))
00255     {
00256       double nSigmaXY       = 100.;
00257       double nSigmaZ        = 100.;
00258       double varFactor      = 4./25.; // Take into account the difference between the RMS and sigma (RMS usually greater than sigma)
00259       double parDistanceXY  = 0.005;  // Unit: [cm]
00260       double parDistanceZ   = 0.5;    // Unit: [cm]
00261       double parDistanceddZ = 1e-3;   // Unit: [rad]
00262       double parDistanceCxy = 1e-5;   // Unit: [cm^2]
00263       double bestEdm        = 1e-1;
00264 
00265       const unsigned int trials = 4;
00266       double largerDist[trials] = {0.1, 5., 10., 100.};
00267 
00268       double covxz,covyz,det;
00269       double deltaMean;
00270       int bestMovementX = 1;
00271       int bestMovementY = 1;
00272       int bestMovementZ = 1;
00273       int goodData;
00274 
00275       double arglist[2];
00276       double amin,errdef,edm;
00277       int nvpar,nparx;
00278       
00279       vector<double>::const_iterator it = vals->begin();
00280 
00281       TFitterMinuit* Gauss3D = new TFitterMinuit(nParams);
00282       if (internalDebug == true) Gauss3D->SetPrintLevel(3);
00283       else Gauss3D->SetPrintLevel(0);
00284       Gauss3D->SetFCN(Gauss3DFunc);
00285       arglist[0] = 10000; // Max number of function calls
00286       arglist[1] = 1e-9;  // Tolerance on likelihood
00287 
00288       if (internalDebug == true) cout << "\n@@@ START FITTING @@@" << endl;
00289 
00290       // @@@ Fit at X-deltaMean | X | X+deltaMean @@@
00291       bestEdm = 1.;
00292       for (int i = 0; i < 3; i++)
00293         {
00294           deltaMean = (double(i)-1.)*std::sqrt((*(it+0))*varFactor);
00295           if (internalDebug == true) cout << "deltaMean --> " << deltaMean << endl;
00296 
00297           Gauss3D->Clear();
00298 
00299           // arg3 - first guess of parameter value
00300           // arg4 - step of the parameter
00301           Gauss3D->SetParameter(0,"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
00302           Gauss3D->SetParameter(1,"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
00303           Gauss3D->SetParameter(2,"var z ", *(it+2), parDistanceZ*parDistanceZ, 0., 0.);
00304           Gauss3D->SetParameter(3,"cov xy", *(it+3), parDistanceCxy, 0., 0.);
00305           Gauss3D->SetParameter(4,"dydz  ", *(it+4), parDistanceddZ, 0., 0.);
00306           Gauss3D->SetParameter(5,"dxdz  ", *(it+5), parDistanceddZ, 0., 0.);
00307           Gauss3D->SetParameter(6,"mean x", *(it+6)+deltaMean, parDistanceXY, 0., 0.);
00308           Gauss3D->SetParameter(7,"mean y", *(it+7), parDistanceXY, 0., 0.);
00309           Gauss3D->SetParameter(8,"mean z", *(it+8), parDistanceZ, 0., 0.);
00310 
00311           // Set the central positions of the centroid for vertex rejection
00312           xPos = Gauss3D->GetParameter(6);
00313           yPos = Gauss3D->GetParameter(7);
00314           zPos = Gauss3D->GetParameter(8);
00315 
00316           // Set dimensions of the centroid for vertex rejection
00317           maxTransRadius = nSigmaXY * std::sqrt(std::fabs(Gauss3D->GetParameter(0)) + std::fabs(Gauss3D->GetParameter(1))) / 2.;
00318           maxLongLength  = nSigmaZ  * std::sqrt(std::fabs(Gauss3D->GetParameter(2)));
00319 
00320           goodData = Gauss3D->ExecuteCommand("MIGRAD",arglist,2);
00321           Gauss3D->GetStats(amin, edm, errdef, nvpar, nparx);
00322 
00323           if (counterVx < minNentries) goodData = -2;
00324           else if (std::isnan(edm) == true) goodData = -1;
00325           else for (unsigned int j = 0; j < nParams; j++) if (std::isnan(Gauss3D->GetParError(j)) == true) { goodData = -1; break; }
00326           if (goodData == 0)
00327             {
00328               covyz = Gauss3D->GetParameter(4)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(1))) - Gauss3D->GetParameter(5)*Gauss3D->GetParameter(3);
00329               covxz = Gauss3D->GetParameter(5)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(0))) - Gauss3D->GetParameter(4)*Gauss3D->GetParameter(3);
00330               
00331               det = std::fabs(Gauss3D->GetParameter(0)) * (std::fabs(Gauss3D->GetParameter(1))*std::fabs(Gauss3D->GetParameter(2)) - covyz*covyz) -
00332                 Gauss3D->GetParameter(3) * (Gauss3D->GetParameter(3)*std::fabs(Gauss3D->GetParameter(2)) - covxz*covyz) +
00333                 covxz * (Gauss3D->GetParameter(3)*covyz - covxz*std::fabs(Gauss3D->GetParameter(1)));
00334               if (det < 0.) { goodData = -1; if (internalDebug == true) cout << "Negative determinant !" << endl; }
00335             }
00336 
00337           if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementX = i; }
00338         }
00339       if (internalDebug == true) cout << "Found bestMovementX --> " << bestMovementX << endl;
00340 
00341       // @@@ Fit at Y-deltaMean | Y | Y+deltaMean @@@
00342       bestEdm = 1.;
00343       for (int i = 0; i < 3; i++)
00344         {
00345           deltaMean = (double(i)-1.)*std::sqrt((*(it+1))*varFactor);
00346           if (internalDebug == true)
00347             {
00348               cout << "deltaMean --> " << deltaMean << endl;
00349               cout << "deltaMean X --> " << (double(bestMovementX)-1.)*std::sqrt((*(it+0))*varFactor) << endl;
00350             }
00351 
00352           Gauss3D->Clear();
00353 
00354           // arg3 - first guess of parameter value
00355           // arg4 - step of the parameter
00356           Gauss3D->SetParameter(0,"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
00357           Gauss3D->SetParameter(1,"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
00358           Gauss3D->SetParameter(2,"var z ", *(it+2), parDistanceZ*parDistanceZ, 0., 0.);
00359           Gauss3D->SetParameter(3,"cov xy", *(it+3), parDistanceCxy, 0., 0.);
00360           Gauss3D->SetParameter(4,"dydz  ", *(it+4), parDistanceddZ, 0., 0.);
00361           Gauss3D->SetParameter(5,"dxdz  ", *(it+5), parDistanceddZ, 0., 0.);
00362           Gauss3D->SetParameter(6,"mean x", *(it+6)+(double(bestMovementX)-1.)*std::sqrt((*(it+0))*varFactor), parDistanceXY, 0., 0.);
00363           Gauss3D->SetParameter(7,"mean y", *(it+7)+deltaMean, parDistanceXY, 0., 0.);
00364           Gauss3D->SetParameter(8,"mean z", *(it+8), parDistanceZ, 0., 0.);
00365 
00366           // Set the central positions of the centroid for vertex rejection
00367           xPos = Gauss3D->GetParameter(6);
00368           yPos = Gauss3D->GetParameter(7);
00369           zPos = Gauss3D->GetParameter(8);
00370 
00371           // Set dimensions of the centroid for vertex rejection
00372           maxTransRadius = nSigmaXY * std::sqrt(std::fabs(Gauss3D->GetParameter(0)) + std::fabs(Gauss3D->GetParameter(1))) / 2.;
00373           maxLongLength  = nSigmaZ  * std::sqrt(std::fabs(Gauss3D->GetParameter(2)));
00374 
00375           goodData = Gauss3D->ExecuteCommand("MIGRAD",arglist,2);
00376           Gauss3D->GetStats(amin, edm, errdef, nvpar, nparx);
00377 
00378           if (counterVx < minNentries) goodData = -2;
00379           else if (std::isnan(edm) == true) goodData = -1;
00380           else for (unsigned int j = 0; j < nParams; j++) if (std::isnan(Gauss3D->GetParError(j)) == true) { goodData = -1; break; }
00381           if (goodData == 0)
00382             {
00383               covyz = Gauss3D->GetParameter(4)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(1))) - Gauss3D->GetParameter(5)*Gauss3D->GetParameter(3);
00384               covxz = Gauss3D->GetParameter(5)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(0))) - Gauss3D->GetParameter(4)*Gauss3D->GetParameter(3);
00385               
00386               det = std::fabs(Gauss3D->GetParameter(0)) * (std::fabs(Gauss3D->GetParameter(1))*std::fabs(Gauss3D->GetParameter(2)) - covyz*covyz) -
00387                 Gauss3D->GetParameter(3) * (Gauss3D->GetParameter(3)*std::fabs(Gauss3D->GetParameter(2)) - covxz*covyz) +
00388                 covxz * (Gauss3D->GetParameter(3)*covyz - covxz*std::fabs(Gauss3D->GetParameter(1)));
00389               if (det < 0.) { goodData = -1; if (internalDebug == true) cout << "Negative determinant !" << endl; }
00390             }
00391           
00392           if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementY = i; }
00393         }
00394       if (internalDebug == true) cout << "Found bestMovementY --> " << bestMovementY << endl;
00395 
00396       // @@@ Fit at Z-deltaMean | Z | Z+deltaMean @@@
00397       bestEdm = 1.;
00398       for (int i = 0; i < 3; i++)
00399         {
00400           deltaMean = (double(i)-1.)*std::sqrt(*(it+2));
00401           if (internalDebug == true)
00402             {
00403               cout << "deltaMean --> " << deltaMean << endl;
00404               cout << "deltaMean X --> " << (double(bestMovementX)-1.)*std::sqrt((*(it+0))*varFactor) << endl;
00405               cout << "deltaMean Y --> " << (double(bestMovementY)-1.)*std::sqrt((*(it+1))*varFactor) << endl;
00406             }
00407 
00408           Gauss3D->Clear();
00409 
00410           // arg3 - first guess of parameter value
00411           // arg4 - step of the parameter
00412           Gauss3D->SetParameter(0,"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
00413           Gauss3D->SetParameter(1,"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
00414           Gauss3D->SetParameter(2,"var z ", *(it+2), parDistanceZ*parDistanceZ, 0., 0.);
00415           Gauss3D->SetParameter(3,"cov xy", *(it+3), parDistanceCxy, 0., 0.);
00416           Gauss3D->SetParameter(4,"dydz  ", *(it+4), parDistanceddZ, 0., 0.);
00417           Gauss3D->SetParameter(5,"dxdz  ", *(it+5), parDistanceddZ, 0., 0.);
00418           Gauss3D->SetParameter(6,"mean x", *(it+6)+(double(bestMovementX)-1.)*std::sqrt((*(it+0))*varFactor), parDistanceXY, 0., 0.);
00419           Gauss3D->SetParameter(7,"mean y", *(it+7)+(double(bestMovementY)-1.)*std::sqrt((*(it+1))*varFactor), parDistanceXY, 0., 0.);
00420           Gauss3D->SetParameter(8,"mean z", *(it+8)+deltaMean, parDistanceZ, 0., 0.);
00421 
00422           // Set the central positions of the centroid for vertex rejection
00423           xPos = Gauss3D->GetParameter(6);
00424           yPos = Gauss3D->GetParameter(7);
00425           zPos = Gauss3D->GetParameter(8);
00426 
00427           // Set dimensions of the centroid for vertex rejection
00428           maxTransRadius = nSigmaXY * std::sqrt(std::fabs(Gauss3D->GetParameter(0)) + std::fabs(Gauss3D->GetParameter(1))) / 2.;
00429           maxLongLength  = nSigmaZ  * std::sqrt(std::fabs(Gauss3D->GetParameter(2)));
00430 
00431           goodData = Gauss3D->ExecuteCommand("MIGRAD",arglist,2);
00432           Gauss3D->GetStats(amin, edm, errdef, nvpar, nparx);
00433 
00434           if (counterVx < minNentries) goodData = -2;
00435           else if (std::isnan(edm) == true) goodData = -1;
00436           else for (unsigned int j = 0; j < nParams; j++) if (std::isnan(Gauss3D->GetParError(j)) == true) { goodData = -1; break; }
00437           if (goodData == 0)
00438             {
00439               covyz = Gauss3D->GetParameter(4)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(1))) - Gauss3D->GetParameter(5)*Gauss3D->GetParameter(3);
00440               covxz = Gauss3D->GetParameter(5)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(0))) - Gauss3D->GetParameter(4)*Gauss3D->GetParameter(3);
00441               
00442               det = std::fabs(Gauss3D->GetParameter(0)) * (std::fabs(Gauss3D->GetParameter(1))*std::fabs(Gauss3D->GetParameter(2)) - covyz*covyz) -
00443                 Gauss3D->GetParameter(3) * (Gauss3D->GetParameter(3)*std::fabs(Gauss3D->GetParameter(2)) - covxz*covyz) +
00444                 covxz * (Gauss3D->GetParameter(3)*covyz - covxz*std::fabs(Gauss3D->GetParameter(1)));
00445               if (det < 0.) { goodData = -1; if (internalDebug == true) cout << "Negative determinant !" << endl; }
00446             }
00447           
00448           if ((goodData == 0) && (std::fabs(edm) < bestEdm)) { bestEdm = edm; bestMovementZ = i; }
00449         }
00450       if (internalDebug == true) cout << "Found bestMovementZ --> " << bestMovementZ << endl;
00451 
00452       Gauss3D->Clear();
00453 
00454       // @@@ FINAL FIT @@@
00455       // arg3 - first guess of parameter value
00456       // arg4 - step of the parameter
00457       Gauss3D->SetParameter(0,"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
00458       Gauss3D->SetParameter(1,"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY, 0., 0.);
00459       Gauss3D->SetParameter(2,"var z ", *(it+2), parDistanceZ*parDistanceZ, 0., 0.);
00460       Gauss3D->SetParameter(3,"cov xy", *(it+3), parDistanceCxy, 0., 0.);
00461       Gauss3D->SetParameter(4,"dydz  ", *(it+4), parDistanceddZ, 0., 0.);
00462       Gauss3D->SetParameter(5,"dxdz  ", *(it+5), parDistanceddZ, 0., 0.);
00463       Gauss3D->SetParameter(6,"mean x", *(it+6)+(double(bestMovementX)-1.)*std::sqrt((*(it+0))*varFactor), parDistanceXY, 0., 0.);
00464       Gauss3D->SetParameter(7,"mean y", *(it+7)+(double(bestMovementY)-1.)*std::sqrt((*(it+1))*varFactor), parDistanceXY, 0., 0.);
00465       Gauss3D->SetParameter(8,"mean z", *(it+8)+(double(bestMovementZ)-1.)*std::sqrt(*(it+2)), parDistanceZ, 0., 0.);
00466 
00467       // Set the central positions of the centroid for vertex rejection
00468       xPos = Gauss3D->GetParameter(6);
00469       yPos = Gauss3D->GetParameter(7);
00470       zPos = Gauss3D->GetParameter(8);
00471       
00472       // Set dimensions of the centroid for vertex rejection
00473       maxTransRadius = nSigmaXY * std::sqrt(std::fabs(Gauss3D->GetParameter(0)) + std::fabs(Gauss3D->GetParameter(1))) / 2.;
00474       maxLongLength  = nSigmaZ  * std::sqrt(std::fabs(Gauss3D->GetParameter(2)));
00475 
00476       goodData = Gauss3D->ExecuteCommand("MIGRAD",arglist,2);      
00477       Gauss3D->GetStats(amin, edm, errdef, nvpar, nparx);
00478       
00479       if (counterVx < minNentries) goodData = -2;
00480       else if (std::isnan(edm) == true) goodData = -1;
00481       else for (unsigned int j = 0; j < nParams; j++) if (std::isnan(Gauss3D->GetParError(j)) == true) { goodData = -1; break; }
00482       if (goodData == 0)
00483         {
00484           covyz = Gauss3D->GetParameter(4)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(1))) - Gauss3D->GetParameter(5)*Gauss3D->GetParameter(3);
00485           covxz = Gauss3D->GetParameter(5)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(0))) - Gauss3D->GetParameter(4)*Gauss3D->GetParameter(3);
00486           
00487           det = std::fabs(Gauss3D->GetParameter(0)) * (std::fabs(Gauss3D->GetParameter(1))*std::fabs(Gauss3D->GetParameter(2)) - covyz*covyz) -
00488             Gauss3D->GetParameter(3) * (Gauss3D->GetParameter(3)*std::fabs(Gauss3D->GetParameter(2)) - covxz*covyz) +
00489             covxz * (Gauss3D->GetParameter(3)*covyz - covxz*std::fabs(Gauss3D->GetParameter(1)));
00490           if (det < 0.) { goodData = -1; if (internalDebug == true) cout << "Negative determinant !" << endl; }
00491         }
00492 
00493       // @@@ FIT WITH DIFFERENT PARAMETER DISTANCES@@@
00494       // arg3 - first guess of parameter value
00495       // arg4 - step of the parameter
00496       for (unsigned int i = 0; i < trials; i++)
00497         {
00498           if ((goodData != 0) && (goodData != -2))
00499             {
00500               Gauss3D->Clear();
00501           
00502               if (internalDebug == true) cout << "FIT WITH DIFFERENT PARAMETER DISTANCES - STEP " << i+1 << endl;      
00503 
00504               Gauss3D->SetParameter(0,"var x ", *(it+0)*varFactor, parDistanceXY*parDistanceXY * largerDist[i], 0, 0);
00505               Gauss3D->SetParameter(1,"var y ", *(it+1)*varFactor, parDistanceXY*parDistanceXY * largerDist[i], 0, 0);
00506               Gauss3D->SetParameter(2,"var z ", *(it+2), parDistanceZ*parDistanceZ * largerDist[i], 0, 0);
00507               Gauss3D->SetParameter(3,"cov xy", *(it+3), parDistanceCxy * largerDist[i], 0, 0);
00508               Gauss3D->SetParameter(4,"dydz  ", *(it+4), parDistanceddZ * largerDist[i], 0, 0);
00509               Gauss3D->SetParameter(5,"dxdz  ", *(it+5), parDistanceddZ * largerDist[i], 0, 0);
00510               Gauss3D->SetParameter(6,"mean x", *(it+6)+(double(bestMovementX)-1.)*std::sqrt((*(it+0))*varFactor), parDistanceXY * largerDist[i], 0, 0);
00511               Gauss3D->SetParameter(7,"mean y", *(it+7)+(double(bestMovementY)-1.)*std::sqrt((*(it+1))*varFactor), parDistanceXY * largerDist[i], 0, 0);
00512               Gauss3D->SetParameter(8,"mean z", *(it+8)+(double(bestMovementZ)-1.)*std::sqrt(*(it+2)), parDistanceZ * largerDist[i], 0, 0);
00513 
00514               // Set the central positions of the centroid for vertex rejection
00515               xPos = Gauss3D->GetParameter(6);
00516               yPos = Gauss3D->GetParameter(7);
00517               zPos = Gauss3D->GetParameter(8);
00518 
00519               // Set dimensions of the centroid for vertex rejection
00520               maxTransRadius = nSigmaXY * std::sqrt(std::fabs(Gauss3D->GetParameter(0)) + std::fabs(Gauss3D->GetParameter(1))) / 2.;
00521               maxLongLength  = nSigmaZ  * std::sqrt(std::fabs(Gauss3D->GetParameter(2)));
00522 
00523               goodData = Gauss3D->ExecuteCommand("MIGRAD",arglist,2);
00524               Gauss3D->GetStats(amin, edm, errdef, nvpar, nparx);
00525       
00526               if (counterVx < minNentries) goodData = -2;
00527               else if (std::isnan(edm) == true) goodData = -1;
00528               else for (unsigned int j = 0; j < nParams; j++) if (std::isnan(Gauss3D->GetParError(j)) == true) { goodData = -1; break; }
00529               if (goodData == 0)
00530                 {
00531                   covyz = Gauss3D->GetParameter(4)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(1))) - Gauss3D->GetParameter(5)*Gauss3D->GetParameter(3);
00532                   covxz = Gauss3D->GetParameter(5)*(std::fabs(Gauss3D->GetParameter(2))-std::fabs(Gauss3D->GetParameter(0))) - Gauss3D->GetParameter(4)*Gauss3D->GetParameter(3);
00533               
00534                   det = std::fabs(Gauss3D->GetParameter(0)) * (std::fabs(Gauss3D->GetParameter(1))*std::fabs(Gauss3D->GetParameter(2)) - covyz*covyz) -
00535                     Gauss3D->GetParameter(3) * (Gauss3D->GetParameter(3)*std::fabs(Gauss3D->GetParameter(2)) - covxz*covyz) +
00536                     covxz * (Gauss3D->GetParameter(3)*covyz - covxz*std::fabs(Gauss3D->GetParameter(1)));
00537                   if (det < 0.) { goodData = -1; if (internalDebug == true) cout << "Negative determinant !" << endl; }
00538                 }
00539             } else break;
00540         }
00541 
00542       if (goodData == 0)
00543         for (unsigned int i = 0; i < nParams; i++)
00544           {
00545             vals->operator[](i) = Gauss3D->GetParameter(i);
00546             vals->operator[](i+nParams) = Gauss3D->GetParError(i);
00547           }
00548       
00549       delete Gauss3D;
00550       return goodData;
00551     }
00552   
00553   return -1;
00554 }
00555 
00556 
00557 void Vx3DHLTAnalyzer::reset(string ResetType)
00558 {
00559   if (ResetType.compare("scratch") == 0)
00560     {
00561       runNumber        = 0;
00562       numberGoodFits   = 0;
00563       numberFits       = 0;
00564       lastLumiOfFit    = 0;
00565       
00566       Vx_X->Reset();
00567       Vx_Y->Reset();
00568       Vx_Z->Reset();
00569       
00570       Vx_ZX->Reset();
00571       Vx_ZY->Reset();
00572       Vx_XY->Reset();
00573       
00574       mXlumi->Reset();
00575       mYlumi->Reset();
00576       mZlumi->Reset();
00577 
00578       sXlumi->Reset();
00579       sYlumi->Reset();
00580       sZlumi->Reset();
00581 
00582       dxdzlumi->Reset();
00583       dydzlumi->Reset();
00584 
00585       hitCounter->Reset();
00586       hitCountHistory->Reset();
00587       goodVxCounter->Reset();
00588       goodVxCountHistory->Reset();
00589       fitResults->Reset();
00590 
00591       reportSummary->Fill(0.);
00592       reportSummaryMap->Fill(0.5, 0.5, 0.);
00593 
00594       Vertices.clear();
00595       
00596       lumiCounter      = 0;
00597       lumiCounterHisto = 0;
00598       totalHits        = 0;
00599       beginTimeOfFit   = 0;
00600       endTimeOfFit     = 0;
00601       beginLumiOfFit   = 0;
00602       endLumiOfFit     = 0;
00603     }
00604   else if (ResetType.compare("whole") == 0)
00605     {
00606       Vx_X->Reset();
00607       Vx_Y->Reset();
00608       Vx_Z->Reset();
00609       
00610       Vx_ZX->Reset();
00611       Vx_ZY->Reset();
00612       Vx_XY->Reset();
00613       
00614       Vertices.clear();
00615       
00616       lumiCounter      = 0;
00617       lumiCounterHisto = 0;
00618       totalHits        = 0;
00619       beginTimeOfFit   = 0;
00620       endTimeOfFit     = 0;
00621       beginLumiOfFit   = 0;
00622       endLumiOfFit     = 0;
00623     }
00624   else if (ResetType.compare("partial") == 0)
00625     {
00626       Vx_X->Reset();
00627       Vx_Y->Reset();
00628       Vx_Z->Reset();
00629       
00630       Vertices.clear();
00631       
00632       lumiCounter      = 0;
00633       totalHits        = 0;
00634       beginTimeOfFit   = 0;
00635       endTimeOfFit     = 0;
00636       beginLumiOfFit   = 0;
00637       endLumiOfFit     = 0;
00638     }
00639   else if (ResetType.compare("nohisto") == 0)
00640     {
00641       Vertices.clear();
00642       
00643       lumiCounter      = 0;
00644       lumiCounterHisto = 0;
00645       totalHits        = 0;
00646       beginTimeOfFit   = 0;
00647       endTimeOfFit     = 0;
00648       beginLumiOfFit   = 0;
00649       endLumiOfFit     = 0;
00650     }
00651   else if (ResetType.compare("hitCounter") == 0)
00652     totalHits          = 0;
00653 }
00654 
00655 
00656 void Vx3DHLTAnalyzer::writeToFile(vector<double>* vals,
00657                                   edm::TimeValue_t BeginTimeOfFit,
00658                                   edm::TimeValue_t EndTimeOfFit,
00659                                   unsigned int BeginLumiOfFit,
00660                                   unsigned int EndLumiOfFit,
00661                                   int dataType)
00662 {
00663   stringstream BufferString;
00664   BufferString.precision(5);
00665 
00666   outputFile.open(fileName.c_str(), ios::out);
00667 
00668   if ((outputFile.is_open() == true) && (vals != NULL) && (vals->size() == 8*2))
00669     {
00670       vector<double>::const_iterator it = vals->begin();
00671 
00672       outputFile << "Runnumber " << runNumber << endl;
00673       outputFile << "BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) << " " << (beginTimeOfFit >> 32) << endl;
00674       outputFile << "EndTimeOfFit " << formatTime(endTimeOfFit >> 32) << " " << (endTimeOfFit >> 32) << endl;
00675       outputFile << "LumiRange " << beginLumiOfFit << " - " << endLumiOfFit << endl;
00676       outputFile << "Type " << dataType << endl;
00677       // 3D Vertexing with Pixel Tracks:
00678       // Good data = Type  3
00679       // Bad data  = Type -1
00680 
00681       BufferString << *(it+0);
00682       outputFile << "X0 " << BufferString.str().c_str() << endl;
00683       BufferString.str("");
00684 
00685       BufferString << *(it+1);
00686       outputFile << "Y0 " << BufferString.str().c_str() << endl;
00687       BufferString.str("");
00688 
00689       BufferString << *(it+2);
00690       outputFile << "Z0 " << BufferString.str().c_str() << endl;
00691       BufferString.str("");
00692 
00693       BufferString << *(it+3);
00694       outputFile << "sigmaZ0 " << BufferString.str().c_str() << endl;
00695       BufferString.str("");
00696 
00697       BufferString << *(it+4);
00698       outputFile << "dxdz " << BufferString.str().c_str() << endl;
00699       BufferString.str("");
00700 
00701       BufferString << *(it+5);
00702       outputFile << "dydz " << BufferString.str().c_str() << endl;
00703       BufferString.str("");
00704 
00705       BufferString << *(it+6);
00706       outputFile << "BeamWidthX " << BufferString.str().c_str() << endl;
00707       BufferString.str("");
00708 
00709       BufferString << *(it+7);
00710       outputFile << "BeamWidthY " << BufferString.str().c_str() << endl;
00711       BufferString.str("");
00712 
00713       outputFile << "Cov(0,j) " << *(it+8) << " 0.0 0.0 0.0 0.0 0.0 0.0" << endl;
00714       outputFile << "Cov(1,j) 0.0 " << *(it+9) << " 0.0 0.0 0.0 0.0 0.0" << endl;
00715       outputFile << "Cov(2,j) 0.0 0.0 " << *(it+10) << " 0.0 0.0 0.0 0.0" << endl;
00716       outputFile << "Cov(3,j) 0.0 0.0 0.0 " << *(it+11) << " 0.0 0.0 0.0" << endl;
00717       outputFile << "Cov(4,j) 0.0 0.0 0.0 0.0 " << *(it+12) << " 0.0 0.0" << endl;
00718       outputFile << "Cov(5,j) 0.0 0.0 0.0 0.0 0.0 " << *(it+13) << " 0.0" << endl;
00719       outputFile << "Cov(6,j) 0.0 0.0 0.0 0.0 0.0 0.0 " << ((*(it+14)) + (*(it+15)) + 2.*std::sqrt((*(it+14))*(*(it+15)))) / 4. << endl;
00720 
00721       outputFile << "EmittanceX 0.0" << endl;
00722       outputFile << "EmittanceY 0.0" << endl;
00723       outputFile << "BetaStar 0.0" << endl;
00724     }
00725   outputFile.close();
00726 
00727   if ((debugMode == true) && (outputDebugFile.is_open() == true) && (vals != NULL) && (vals->size() == 8*2))
00728     {
00729       vector<double>::const_iterator it = vals->begin();
00730 
00731       outputDebugFile << "Runnumber " << runNumber << endl;
00732       outputDebugFile << "BeginTimeOfFit " << formatTime(beginTimeOfFit >> 32) << " " << (beginTimeOfFit >> 32) << endl;
00733       outputDebugFile << "EndTimeOfFit " << formatTime(endTimeOfFit >> 32) << " " << (endTimeOfFit >> 32) << endl;
00734       outputDebugFile << "LumiRange " << beginLumiOfFit << " - " << endLumiOfFit << endl;
00735       outputDebugFile << "Type " << dataType << endl;
00736       // 3D Vertexing with Pixel Tracks:
00737       // Good data = Type  3
00738       // Bad data  = Type -1
00739           
00740       BufferString << *(it+0);
00741       outputDebugFile << "X0 " << BufferString.str().c_str() << endl;
00742       BufferString.str("");
00743           
00744       BufferString << *(it+1);
00745       outputDebugFile << "Y0 " << BufferString.str().c_str() << endl;
00746       BufferString.str("");
00747           
00748       BufferString << *(it+2);
00749       outputDebugFile << "Z0 " << BufferString.str().c_str() << endl;
00750       BufferString.str("");
00751           
00752       BufferString << *(it+3);
00753       outputDebugFile << "sigmaZ0 " << BufferString.str().c_str() << endl;
00754       BufferString.str("");
00755           
00756       BufferString << *(it+4);
00757       outputDebugFile << "dxdz " << BufferString.str().c_str() << endl;
00758       BufferString.str("");
00759           
00760       BufferString << *(it+5);
00761       outputDebugFile << "dydz " << BufferString.str().c_str() << endl;
00762       BufferString.str("");
00763           
00764       BufferString << *(it+6);
00765       outputDebugFile << "BeamWidthX " << BufferString.str().c_str() << endl;
00766       BufferString.str("");
00767           
00768       BufferString << *(it+7);
00769       outputDebugFile << "BeamWidthY " << BufferString.str().c_str() << endl;
00770       BufferString.str("");
00771           
00772       outputDebugFile << "Cov(0,j) " << *(it+8) << " 0.0 0.0 0.0 0.0 0.0 0.0" << endl;
00773       outputDebugFile << "Cov(1,j) 0.0 " << *(it+9) << " 0.0 0.0 0.0 0.0 0.0" << endl;
00774       outputDebugFile << "Cov(2,j) 0.0 0.0 " << *(it+10) << " 0.0 0.0 0.0 0.0" << endl;
00775       outputDebugFile << "Cov(3,j) 0.0 0.0 0.0 " << *(it+11) << " 0.0 0.0 0.0" << endl;
00776       outputDebugFile << "Cov(4,j) 0.0 0.0 0.0 0.0 " << *(it+12) << " 0.0 0.0" << endl;
00777       outputDebugFile << "Cov(5,j) 0.0 0.0 0.0 0.0 0.0 " << *(it+13) << " 0.0" << endl;
00778       outputDebugFile << "Cov(6,j) 0.0 0.0 0.0 0.0 0.0 0.0 " << ((*(it+14)) + (*(it+15)) + 2.*std::sqrt((*(it+14))*(*(it+15)))) / 4. << endl;
00779           
00780       outputDebugFile << "EmittanceX 0.0" << endl;
00781       outputDebugFile << "EmittanceY 0.0" << endl;
00782       outputDebugFile << "BetaStar 0.0" << endl;
00783     }
00784 }
00785 
00786 
00787 void Vx3DHLTAnalyzer::beginLuminosityBlock(const LuminosityBlock& lumiBlock, 
00788                                            const EventSetup& iSetup)
00789 {
00790   if ((lumiCounter == 0) && (lumiBlock.luminosityBlock() > lastLumiOfFit))
00791     {
00792       beginTimeOfFit = lumiBlock.beginTime().value();
00793       beginLumiOfFit = lumiBlock.luminosityBlock();
00794       lumiCounter++;
00795       lumiCounterHisto++;
00796     }
00797   else if ((lumiCounter != 0) && (lumiBlock.luminosityBlock() >= (beginLumiOfFit+lumiCounter))) { lumiCounter++; lumiCounterHisto++; }
00798 }
00799 
00800 
00801 void Vx3DHLTAnalyzer::endLuminosityBlock(const LuminosityBlock& lumiBlock,
00802                                          const EventSetup& iSetup)
00803 {
00804   stringstream histTitle;
00805   int goodData;
00806   unsigned int nParams = 9;
00807 
00808   if ((lumiCounter%nLumiReset == 0) && (nLumiReset != 0) && (beginTimeOfFit != 0) && (runNumber != 0))
00809     {
00810       endTimeOfFit  = lumiBlock.endTime().value();
00811       endLumiOfFit  = lumiBlock.luminosityBlock();
00812       lastLumiOfFit = endLumiOfFit;
00813       vector<double> vals;
00814 
00815       hitCounter->ShiftFillLast((double)totalHits, std::sqrt((double)totalHits), nLumiReset);
00816 
00817       if (lastLumiOfFit % prescaleHistory == 0)
00818         {
00819           hitCountHistory->getTH1()->SetBinContent(lastLumiOfFit, (double)totalHits);
00820           hitCountHistory->getTH1()->SetBinError(lastLumiOfFit, std::sqrt((double)totalHits));
00821         }
00822 
00823       if (dataFromFit == true)
00824         {
00825           vector<double> fitResults;
00826 
00827           fitResults.push_back(Vx_X->getTH1()->GetRMS()*Vx_X->getTH1()->GetRMS());
00828           fitResults.push_back(Vx_Y->getTH1()->GetRMS()*Vx_Y->getTH1()->GetRMS());
00829           fitResults.push_back(Vx_Z->getTH1()->GetRMS()*Vx_Z->getTH1()->GetRMS());
00830           fitResults.push_back(0.0);
00831           fitResults.push_back(0.0);
00832           fitResults.push_back(0.0);
00833           fitResults.push_back(Vx_X->getTH1()->GetMean());
00834           fitResults.push_back(Vx_Y->getTH1()->GetMean());
00835           fitResults.push_back(Vx_Z->getTH1()->GetMean());
00836           for (unsigned int i = 0; i < nParams; i++) fitResults.push_back(0.0);
00837           
00838           goodData = MyFit(&fitResults);                      
00839 
00840           if (internalDebug == true) 
00841             {
00842               cout << "goodData --> " << goodData << endl;
00843               cout << "Used vertices --> " << counterVx << endl;
00844               cout << "var x -->  " << fitResults[0] << " +/- " << fitResults[0+nParams] << endl;
00845               cout << "var y -->  " << fitResults[1] << " +/- " << fitResults[1+nParams] << endl;
00846               cout << "var z -->  " << fitResults[2] << " +/- " << fitResults[2+nParams] << endl;
00847               cout << "cov xy --> " << fitResults[3] << " +/- " << fitResults[3+nParams] << endl;
00848               cout << "dydz   --> " << fitResults[4] << " +/- " << fitResults[4+nParams] << endl;
00849               cout << "dxdz   --> " << fitResults[5] << " +/- " << fitResults[5+nParams] << endl;
00850               cout << "mean x --> " << fitResults[6] << " +/- " << fitResults[6+nParams] << endl;
00851               cout << "mean y --> " << fitResults[7] << " +/- " << fitResults[7+nParams] << endl;
00852               cout << "mean z --> " << fitResults[8] << " +/- " << fitResults[8+nParams] << endl;
00853             }
00854 
00855           if (goodData == 0)
00856             {            
00857               vals.push_back(fitResults[6]);
00858               vals.push_back(fitResults[7]);
00859               vals.push_back(fitResults[8]);
00860               vals.push_back(std::sqrt(std::fabs(fitResults[2])));
00861               vals.push_back(fitResults[5]);
00862               vals.push_back(fitResults[4]);
00863               vals.push_back(std::sqrt(std::fabs(fitResults[0])));
00864               vals.push_back(std::sqrt(std::fabs(fitResults[1])));
00865 
00866               vals.push_back(powf(fitResults[6+nParams],2.));
00867               vals.push_back(powf(fitResults[7+nParams],2.));
00868               vals.push_back(powf(fitResults[8+nParams],2.));
00869               vals.push_back(powf(std::fabs(fitResults[2+nParams]) / (2.*std::sqrt(std::fabs(fitResults[2]))),2.));
00870               vals.push_back(powf(fitResults[5+nParams],2.));
00871               vals.push_back(powf(fitResults[4+nParams],2.));
00872               vals.push_back(powf(std::fabs(fitResults[0+nParams]) / (2.*std::sqrt(std::fabs(fitResults[0]))),2.));
00873               vals.push_back(powf(std::fabs(fitResults[1+nParams]) / (2.*std::sqrt(std::fabs(fitResults[1]))),2.));
00874             }
00875           else for (unsigned int i = 0; i < 8*2; i++) vals.push_back(0.0);
00876 
00877           fitResults.clear();
00878         }
00879       else
00880         {
00881           counterVx = Vx_X->getTH1F()->GetEntries();
00882             
00883           if (Vx_X->getTH1F()->GetEntries() >= minNentries)
00884             {
00885             goodData = 0;
00886             
00887             vals.push_back(Vx_X->getTH1F()->GetMean());
00888             vals.push_back(Vx_Y->getTH1F()->GetMean());
00889             vals.push_back(Vx_Z->getTH1F()->GetMean());
00890             vals.push_back(Vx_Z->getTH1F()->GetRMS());
00891             vals.push_back(0.0);
00892             vals.push_back(0.0);
00893             vals.push_back(Vx_X->getTH1F()->GetRMS());
00894             vals.push_back(Vx_Y->getTH1F()->GetRMS());
00895             
00896             vals.push_back(powf(Vx_X->getTH1F()->GetMeanError(),2.));
00897             vals.push_back(powf(Vx_Y->getTH1F()->GetMeanError(),2.));
00898             vals.push_back(powf(Vx_Z->getTH1F()->GetMeanError(),2.));
00899             vals.push_back(powf(Vx_Z->getTH1F()->GetRMSError(),2.));
00900             vals.push_back(0.0);
00901             vals.push_back(0.0);
00902             vals.push_back(powf(Vx_X->getTH1F()->GetRMSError(),2.));
00903             vals.push_back(powf(Vx_Y->getTH1F()->GetRMSError(),2.));
00904             }
00905           else
00906             {
00907               goodData = -2;
00908               for (unsigned int i = 0; i < 8*2; i++) vals.push_back(0.0);
00909             }
00910         }
00911 
00912       // vals[0]  = X0
00913       // vals[1]  = Y0
00914       // vals[2]  = Z0
00915       // vals[3]  = sigmaZ0
00916       // vals[4]  = dxdz
00917       // vals[5]  = dydz
00918       // vals[6]  = BeamWidthX
00919       // vals[7]  = BeamWidthY
00920 
00921       // vals[8]  = err^2 X0
00922       // vals[9]  = err^2 Y0
00923       // vals[10] = err^2 Z0
00924       // vals[11] = err^2 sigmaZ0
00925       // vals[12] = err^2 dxdz
00926       // vals[13] = err^2 dydz
00927       // vals[14] = err^2 BeamWidthX
00928       // vals[15] = err^2 BeamWidthY
00929 
00930       // "goodData" CODE:
00931       //  0 == OK --> Reset
00932       // -2 == NO OK - not enough "minNentries" --> Wait for more lumisections
00933       // Any other number == NO OK --> Reset
00934 
00935       numberFits++;
00936       if (goodData == 0)
00937         {
00938           writeToFile(&vals, beginTimeOfFit, endTimeOfFit, beginLumiOfFit, endLumiOfFit, 3);
00939           if ((internalDebug == true) && (outputDebugFile.is_open() == true)) outputDebugFile << "Used vertices: " << counterVx << endl;
00940 
00941           numberGoodFits++;
00942 
00943           histTitle << "Fitted Beam Spot [cm] (Lumi start: " << beginLumiOfFit << " - Lumi end: " << endLumiOfFit << ")";
00944           if (lumiCounterHisto >= maxLumiIntegration) reset("whole");
00945           else reset("partial");
00946         }
00947       else
00948         {
00949           writeToFile(&vals, beginTimeOfFit, endTimeOfFit, beginLumiOfFit, endLumiOfFit, -1);
00950           if ((internalDebug == true) && (outputDebugFile.is_open() == true)) outputDebugFile << "Used vertices: " << counterVx << endl;
00951 
00952           if (goodData == -2)
00953             {
00954               histTitle << "Fitted Beam Spot [cm] (not enough statistics)";
00955               if (lumiCounter >= maxLumiIntegration) reset("whole");
00956               else reset("hitCounter");
00957             }
00958           else
00959             {
00960               histTitle << "Fitted Beam Spot [cm] (problems)";
00961               if (lumiCounterHisto >= maxLumiIntegration) reset("whole");
00962               else reset("partial");
00963 
00964               counterVx = 0;
00965             }
00966         }
00967 
00968       reportSummary->Fill(numberFits != 0 ? (double)numberGoodFits/(double)numberFits : 0.0);
00969       reportSummaryMap->Fill(0.5, 0.5, numberFits != 0 ? (double)numberGoodFits/(double)numberFits : 0.0);
00970 
00971       fitResults->setAxisTitle(histTitle.str().c_str(), 1);
00972       
00973       fitResults->setBinContent(1, 9, vals[0]);
00974       fitResults->setBinContent(1, 8, vals[1]);
00975       fitResults->setBinContent(1, 7, vals[2]);
00976       fitResults->setBinContent(1, 6, vals[3]);
00977       fitResults->setBinContent(1, 5, vals[4]);
00978       fitResults->setBinContent(1, 4, vals[5]);
00979       fitResults->setBinContent(1, 3, vals[6]);
00980       fitResults->setBinContent(1, 2, vals[7]);
00981       fitResults->setBinContent(1, 1, counterVx);
00982       
00983       fitResults->setBinContent(2, 9, std::sqrt(vals[8]));
00984       fitResults->setBinContent(2, 8, std::sqrt(vals[9]));
00985       fitResults->setBinContent(2, 7, std::sqrt(vals[10]));
00986       fitResults->setBinContent(2, 6, std::sqrt(vals[11]));
00987       fitResults->setBinContent(2, 5, std::sqrt(vals[12]));
00988       fitResults->setBinContent(2, 4, std::sqrt(vals[13]));
00989       fitResults->setBinContent(2, 3, std::sqrt(vals[14]));
00990       fitResults->setBinContent(2, 2, std::sqrt(vals[15]));
00991       fitResults->setBinContent(2, 1, std::sqrt(counterVx));
00992 
00993       // Linear fit to the historical plots
00994       TF1* myLinFit = new TF1("myLinFit", "[0] + [1]*x", mXlumi->getTH1()->GetXaxis()->GetXmin(), mXlumi->getTH1()->GetXaxis()->GetXmax());
00995       myLinFit->SetLineColor(2);
00996       myLinFit->SetLineWidth(2);
00997       myLinFit->SetParName(0,"Intercept");
00998       myLinFit->SetParName(1,"Slope");
00999 
01000       mXlumi->ShiftFillLast(vals[0], std::sqrt(vals[8]), nLumiReset);
01001       myLinFit->SetParameter(0, mXlumi->getTH1()->GetMean(2));
01002       myLinFit->SetParameter(1, 0.0);
01003       mXlumi->getTH1()->Fit("myLinFit","QR");
01004 
01005       mYlumi->ShiftFillLast(vals[1], std::sqrt(vals[9]), nLumiReset);
01006       myLinFit->SetParameter(0, mYlumi->getTH1()->GetMean(2));
01007       myLinFit->SetParameter(1, 0.0);
01008       mYlumi->getTH1()->Fit("myLinFit","QR");
01009 
01010       mZlumi->ShiftFillLast(vals[2], std::sqrt(vals[10]), nLumiReset);
01011       myLinFit->SetParameter(0, mZlumi->getTH1()->GetMean(2));
01012       myLinFit->SetParameter(1, 0.0);
01013       mZlumi->getTH1()->Fit("myLinFit","QR");
01014 
01015       sXlumi->ShiftFillLast(vals[6], std::sqrt(vals[14]), nLumiReset);
01016       myLinFit->SetParameter(0, sXlumi->getTH1()->GetMean(2));
01017       myLinFit->SetParameter(1, 0.0);
01018       sXlumi->getTH1()->Fit("myLinFit","QR");
01019 
01020       sYlumi->ShiftFillLast(vals[7], std::sqrt(vals[15]), nLumiReset);
01021       myLinFit->SetParameter(0, sYlumi->getTH1()->GetMean(2));
01022       myLinFit->SetParameter(1, 0.0);
01023       sYlumi->getTH1()->Fit("myLinFit","QR");
01024 
01025       sZlumi->ShiftFillLast(vals[3], std::sqrt(vals[11]), nLumiReset);
01026       myLinFit->SetParameter(0, sZlumi->getTH1()->GetMean(2));
01027       myLinFit->SetParameter(1, 0.0);
01028       sZlumi->getTH1()->Fit("myLinFit","QR");
01029 
01030       dxdzlumi->ShiftFillLast(vals[4], std::sqrt(vals[12]), nLumiReset);
01031       myLinFit->SetParameter(0, dxdzlumi->getTH1()->GetMean(2));
01032       myLinFit->SetParameter(1, 0.0);
01033       dxdzlumi->getTH1()->Fit("myLinFit","QR");
01034 
01035       dydzlumi->ShiftFillLast(vals[5], std::sqrt(vals[13]), nLumiReset);
01036       myLinFit->SetParameter(0, dydzlumi->getTH1()->GetMean(2));
01037       myLinFit->SetParameter(1, 0.0);
01038       dydzlumi->getTH1()->Fit("myLinFit","QR");
01039       
01040       goodVxCounter->ShiftFillLast((double)counterVx, std::sqrt((double)counterVx), nLumiReset);      
01041       myLinFit->SetParameter(0, goodVxCounter->getTH1()->GetMean(2));
01042       myLinFit->SetParameter(1, 0.0);
01043       goodVxCounter->getTH1()->Fit("myLinFit","QR");
01044 
01045       if (lastLumiOfFit % prescaleHistory == 0)
01046         {
01047           goodVxCountHistory->getTH1()->SetBinContent(lastLumiOfFit, (double)counterVx);
01048           goodVxCountHistory->getTH1()->SetBinError(lastLumiOfFit, std::sqrt((double)counterVx));
01049         }
01050 
01051       delete myLinFit;
01052 
01053       vals.clear();
01054     }
01055   else if (nLumiReset == 0)
01056     {
01057       histTitle << "Fitted Beam Spot [cm] (no ongoing fits)";
01058       fitResults->setAxisTitle(histTitle.str().c_str(), 1);
01059       reportSummaryMap->Fill(0.5, 0.5, 1.0);
01060       hitCounter->ShiftFillLast(totalHits, std::sqrt(totalHits), 1);
01061       reset("nohisto");
01062     }
01063 }
01064 
01065 
01066 void Vx3DHLTAnalyzer::beginJob()
01067 {
01068   DQMStore* dbe = 0;
01069   dbe = Service<DQMStore>().operator->();
01070  
01071   // ### Set internal variables ###
01072   nBinsHistoricalPlot = 80;
01073   nBinsWholeHistory   = 3000; // Corresponds to about 20h of data taking: 20h * 60min * 60s / 23s per lumi-block = 3130
01074   // ##############################
01075 
01076   if ( dbe ) 
01077     {
01078       dbe->setCurrentFolder("BeamPixel");
01079 
01080       Vx_X = dbe->book1D("vertex x", "Primary Vertex X Coordinate Distribution", int(rint(xRange/xStep)), -xRange/2., xRange/2.);
01081       Vx_Y = dbe->book1D("vertex y", "Primary Vertex Y Coordinate Distribution", int(rint(yRange/yStep)), -yRange/2., yRange/2.);
01082       Vx_Z = dbe->book1D("vertex z", "Primary Vertex Z Coordinate Distribution", int(rint(zRange/zStep)), -zRange/2., zRange/2.);
01083       Vx_X->setAxisTitle("Primary Vertices X [cm]",1);
01084       Vx_X->setAxisTitle("Entries [#]",2);
01085       Vx_Y->setAxisTitle("Primary Vertices Y [cm]",1);
01086       Vx_Y->setAxisTitle("Entries [#]",2);
01087       Vx_Z->setAxisTitle("Primary Vertices Z [cm]",1);
01088       Vx_Z->setAxisTitle("Entries [#]",2);
01089  
01090       mXlumi = dbe->book1D("muX vs lumi", "\\mu_{x} vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5);
01091       mYlumi = dbe->book1D("muY vs lumi", "\\mu_{y} vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5);
01092       mZlumi = dbe->book1D("muZ vs lumi", "\\mu_{z} vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5);
01093       mXlumi->setAxisTitle("Lumisection [#]",1);
01094       mXlumi->setAxisTitle("\\mu_{x} [cm]",2);
01095       mXlumi->getTH1()->SetOption("E1");
01096       mYlumi->setAxisTitle("Lumisection [#]",1);
01097       mYlumi->setAxisTitle("\\mu_{y} [cm]",2);
01098       mYlumi->getTH1()->SetOption("E1");
01099       mZlumi->setAxisTitle("Lumisection [#]",1);
01100       mZlumi->setAxisTitle("\\mu_{z} [cm]",2);
01101       mZlumi->getTH1()->SetOption("E1");
01102 
01103       sXlumi = dbe->book1D("sigmaX vs lumi", "\\sigma_{x} vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5);
01104       sYlumi = dbe->book1D("sigmaY vs lumi", "\\sigma_{y} vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5);
01105       sZlumi = dbe->book1D("sigmaZ vs lumi", "\\sigma_{z} vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5);
01106       sXlumi->setAxisTitle("Lumisection [#]",1);
01107       sXlumi->setAxisTitle("\\sigma_{x} [cm]",2);
01108       sXlumi->getTH1()->SetOption("E1");
01109       sYlumi->setAxisTitle("Lumisection [#]",1);
01110       sYlumi->setAxisTitle("\\sigma_{y} [cm]",2);
01111       sYlumi->getTH1()->SetOption("E1");
01112       sZlumi->setAxisTitle("Lumisection [#]",1);
01113       sZlumi->setAxisTitle("\\sigma_{z} [cm]",2);
01114       sZlumi->getTH1()->SetOption("E1");
01115 
01116       dxdzlumi = dbe->book1D("dxdz vs lumi", "dX/dZ vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5);
01117       dydzlumi = dbe->book1D("dydz vs lumi", "dY/dZ vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5);
01118       dxdzlumi->setAxisTitle("Lumisection [#]",1);
01119       dxdzlumi->setAxisTitle("dX/dZ [rad]",2);
01120       dxdzlumi->getTH1()->SetOption("E1");
01121       dydzlumi->setAxisTitle("Lumisection [#]",1);
01122       dydzlumi->setAxisTitle("dY/dZ [rad]",2);
01123       dydzlumi->getTH1()->SetOption("E1");
01124 
01125       Vx_ZX = dbe->book2D("vertex zx", "Primary Vertex ZX Coordinate Distribution", int(rint(zRange/zStep/5.)), -zRange/2., zRange/2., int(rint(xRange/xStep/5.)), -xRange/2., xRange/2.);
01126       Vx_ZY = dbe->book2D("vertex zy", "Primary Vertex ZY Coordinate Distribution", int(rint(zRange/zStep/5.)), -zRange/2., zRange/2., int(rint(yRange/yStep/5.)), -yRange/2., yRange/2.);
01127       Vx_XY = dbe->book2D("vertex xy", "Primary Vertex XY Coordinate Distribution", int(rint(xRange/xStep/5.)), -xRange/2., xRange/2., int(rint(yRange/yStep/5.)), -yRange/2., yRange/2.);
01128       Vx_ZX->setAxisTitle("Primary Vertices Z [cm]",1);
01129       Vx_ZX->setAxisTitle("Primary Vertices X [cm]",2);
01130       Vx_ZX->setAxisTitle("Entries [#]",3);
01131       Vx_ZY->setAxisTitle("Primary Vertices Z [cm]",1);
01132       Vx_ZY->setAxisTitle("Primary Vertices Y [cm]",2);
01133       Vx_ZY->setAxisTitle("Entries [#]",3);
01134       Vx_XY->setAxisTitle("Primary Vertices X [cm]",1);
01135       Vx_XY->setAxisTitle("Primary Vertices Y [cm]",2);
01136       Vx_XY->setAxisTitle("Entries [#]",3);
01137 
01138       hitCounter = dbe->book1D("pixelHits vs lumi", "# Pixel-Hits vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5);
01139       hitCounter->setAxisTitle("Lumisection [#]",1);
01140       hitCounter->setAxisTitle("Pixel-Hits [#]",2);
01141       hitCounter->getTH1()->SetOption("E1");
01142 
01143       hitCountHistory = dbe->book1D("hist pixelHits vs lumi", "History: # Pixel-Hits vs. Lumi", nBinsWholeHistory, 0.5, (double)nBinsWholeHistory+0.5);
01144       hitCountHistory->setAxisTitle("Lumisection [#]",1);
01145       hitCountHistory->setAxisTitle("Pixel-Hits [#]",2);
01146       hitCountHistory->getTH1()->SetOption("E1");
01147 
01148       goodVxCounter = dbe->book1D("good vertices vs lumi", "# Good vertices vs. Lumisection", nBinsHistoricalPlot, 0.5, (double)nBinsHistoricalPlot+0.5);
01149       goodVxCounter->setAxisTitle("Lumisection [#]",1);
01150       goodVxCounter->setAxisTitle("Good vertices [#]",2);
01151       goodVxCounter->getTH1()->SetOption("E1");
01152 
01153       goodVxCountHistory = dbe->book1D("hist good vx vs lumi", "History: # Good vx vs. Lumi", nBinsWholeHistory, 0.5, (double)nBinsWholeHistory+0.5);
01154       goodVxCountHistory->setAxisTitle("Lumisection [#]",1);
01155       goodVxCountHistory->setAxisTitle("Good vertices [#]",2);
01156       goodVxCountHistory->getTH1()->SetOption("E1");
01157 
01158       fitResults = dbe->book2D("fit results","Results of Beam Spot Fit", 2, 0., 2., 9, 0., 9.);
01159       fitResults->setAxisTitle("Fitted Beam Spot [cm]", 1);
01160       fitResults->setBinLabel(9, "X", 2);
01161       fitResults->setBinLabel(8, "Y", 2);
01162       fitResults->setBinLabel(7, "Z", 2);
01163       fitResults->setBinLabel(6, "\\sigma_{Z}", 2);
01164       fitResults->setBinLabel(5, "#frac{dX}{dZ}[rad]", 2);
01165       fitResults->setBinLabel(4, "#frac{dY}{dZ}[rad]", 2);
01166       fitResults->setBinLabel(3, "\\sigma_{X}", 2);
01167       fitResults->setBinLabel(2, "\\sigma_{Y}", 2);
01168       fitResults->setBinLabel(1, "Vertices", 2);
01169       fitResults->setBinLabel(1, "Value", 1);
01170       fitResults->setBinLabel(2, "Stat. Error", 1);
01171       fitResults->getTH1()->SetOption("text");
01172 
01173       dbe->setCurrentFolder("BeamPixel/EventInfo");
01174       reportSummary = dbe->bookFloat("reportSummary");
01175       reportSummary->Fill(0.);
01176       reportSummaryMap = dbe->book2D("reportSummaryMap","Pixel-Vertices Beam Spot: % Good Fits", 1, 0., 1., 1, 0., 1.);
01177       reportSummaryMap->Fill(0.5, 0.5, 0.);
01178       dbe->setCurrentFolder("BeamPixel/EventInfo/reportSummaryContents");
01179 
01180       // Convention for reportSummary and reportSummaryMap:
01181       // - 0%  at the moment of creation of the histogram
01182       // - n%  numberGoodFits / numberFits
01183     }
01184 
01185   // ### Set internal variables ###
01186   reset("scratch");
01187   prescaleHistory      = 1;
01188   maxLumiIntegration   = 15;
01189   minVxDoF             = 4.;
01190   // For vertex fitter without track-weight: d.o.f. = 2*NTracks - 3
01191   // For vertex fitter with track-weight:    d.o.f. = sum_NTracks(2*track_weight) - 3
01192   internalDebug        = false;
01193   considerVxCovariance = true;
01194   pi = 3.141592653589793238;
01195   // ##############################
01196 }
01197 
01198 
01199 void Vx3DHLTAnalyzer::endJob() { reset("scratch"); }
01200 
01201 
01202 // Define this as a plug-in
01203 DEFINE_FWK_MODULE(Vx3DHLTAnalyzer);
/data/refman/pasoursint/CMSSW_5_3_8_patch3/src/DQM/BeamMonitor/plugins/Vx3DHLTAnalyzer.cc
