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