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