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