CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/DQM/DTMonitorClient/src/DTOccupancyCluster.cc

Go to the documentation of this file.
00001 
00002 /*
00003  *  See header file for a description of this class.
00004  *
00005  *  $Date: 2008/10/16 09:33:39 $
00006  *  $Revision: 1.3 $
00007  *  \author G. Cerminara - INFN Torino
00008  */
00009 
00010 #include "DTOccupancyCluster.h"
00011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00012 
00013 #include "TH2F.h"
00014 #include "TMath.h"
00015 
00016 #include <iostream>
00017 
00018 using namespace std;
00019 using namespace edm;
00020 
00021 
00022 DTOccupancyCluster::DTOccupancyCluster(const DTOccupancyPoint& firstPoint,
00023                                        const DTOccupancyPoint& secondPoint) : radius(0),
00024                                                                               theMaxMean(-1.),
00025                                                                               theMaxRMS(-1.),
00026                                                                               theMeanSum(0.),
00027                                                                               theRMSSum(0.) {
00028   if(!qualityCriterion(firstPoint,secondPoint)) {
00029     theValidity = false;
00030   } else {
00031     // compute the cluster quantities
00032     thePoints.push_back(firstPoint);
00033     thePoints.push_back(secondPoint);
00034     if(firstPoint.mean() > secondPoint.mean()) theMaxMean = firstPoint.mean();
00035     else theMaxMean = secondPoint.mean();
00036 
00037     if(firstPoint.rms() > secondPoint.rms()) theMaxRMS = firstPoint.rms();
00038     else theMaxRMS = secondPoint.rms();
00039     theMeanSum += firstPoint.mean();
00040     theRMSSum += firstPoint.rms();
00041 
00042     theMeanSum += secondPoint.mean();
00043     theRMSSum += secondPoint.rms();
00044 
00045 
00046     computeRadius();
00047   }
00048 }
00049 
00050 
00051 
00052 DTOccupancyCluster::DTOccupancyCluster(const DTOccupancyPoint& singlePoint) : radius(0),
00053                                                                               theMaxMean(singlePoint.mean()),
00054                                                                               theMaxRMS(singlePoint.rms()),
00055                                                                               theMeanSum(singlePoint.mean()),
00056                                                                               theRMSSum(singlePoint.rms()) {
00057   theValidity = true;
00058 
00059   // compute the cluster quantities
00060   thePoints.push_back(singlePoint);
00061 }
00062 
00063 
00064 
00065 DTOccupancyCluster::~DTOccupancyCluster(){}
00066 
00067 
00068 
00069  // Check if the cluster candidate satisfies the quality requirements
00070 bool DTOccupancyCluster::isValid() const {
00071   return theValidity;
00072 }
00073 
00074 
00075 
00076 // Add a point to the cluster: returns false if the point does not satisfy the
00077 // quality requirement
00078 bool DTOccupancyCluster::addPoint(const DTOccupancyPoint& anotherPoint) {
00079   LogTrace("DTDQM|DTMonitorClient|DTOccupancyTest|DTOccupancyCluster")
00080     << "   Add a point to the cluster: mean: " << anotherPoint.mean()
00081     << " rms: " << anotherPoint.rms() << endl;
00082   if(qualityCriterion(anotherPoint)) {
00083     LogTrace("DTDQM|DTMonitorClient|DTOccupancyTest|DTOccupancyCluster") << "   point is valid" << endl;
00084     thePoints.push_back(anotherPoint);
00085     // Compute the new cluster size
00086     computeRadius();
00087     // compute the max mean and RMS
00088     if(anotherPoint.mean() > theMaxMean) {
00089       theMaxMean = anotherPoint.mean();
00090     }
00091     if(anotherPoint.rms() > theMaxRMS) {
00092       theMaxRMS = anotherPoint.rms();
00093     }
00094     theMeanSum += anotherPoint.mean();
00095     theRMSSum += anotherPoint.rms();
00096     return true;
00097   } 
00098   return false;
00099 }
00100 
00101 
00102 
00103 // Compute the distance of a single point from the cluster
00104 // (minimum distance with respect to the cluster points)
00105 double DTOccupancyCluster::distance(const DTOccupancyPoint& point) const {
00106   double dist = 99999999;
00107   // compute the minimum distance from a point
00108   for(vector<DTOccupancyPoint>::const_iterator pt = thePoints.begin();
00109       pt != thePoints.end(); ++pt) {
00110     double distance = point.distance(*pt);
00111     if(distance < dist) {
00112       dist = distance;
00113     }
00114   }
00115   return dist;
00116 }
00117 
00118 
00119 
00120 double  DTOccupancyCluster::averageMean() const {
00121   return theMeanSum/(double)thePoints.size();
00122 }
00123 
00124 
00125 
00126 double  DTOccupancyCluster::averageRMS() const {
00127   return theRMSSum/(double)thePoints.size();
00128 }
00129 
00130 
00131 
00132 double  DTOccupancyCluster::maxMean() const {
00133   return theMaxMean;
00134 }
00135 
00136 
00137   
00138 double  DTOccupancyCluster::maxRMS() const {
00139   return theMaxRMS;
00140 }
00141 
00142 
00143 
00144 TH2F * DTOccupancyCluster::getHisto(std::string histoName, int nBinsX, double minX, double maxX,
00145                                     int nBinsY, double minY, double maxY, int fillColor) const {
00146   TH2F *histo = new TH2F(histoName.c_str(),histoName.c_str(),
00147                          nBinsX, minX, maxX, nBinsY, minY, maxY);
00148   histo->SetFillColor(fillColor);
00149   for(vector<DTOccupancyPoint>::const_iterator pt = thePoints.begin();
00150       pt != thePoints.end(); ++pt) {
00151     histo->Fill((*pt).mean(), (*pt).rms());
00152   }
00153   return histo;
00154 }
00155 
00156 
00157 
00158 bool DTOccupancyCluster::qualityCriterion(const DTOccupancyPoint& firstPoint,
00159                                           const DTOccupancyPoint& secondPoint) {
00160   //   cout << "   Seed qualityCriterion results: " << endl;
00161   //   cout << "      delta mean: " << firstPoint.deltaMean(secondPoint) << endl;
00162   //   cout << "      delta rms: " << firstPoint.deltaRMS(secondPoint) << endl;
00163   //   cout << "      average rms: " << computeAverageRMS(firstPoint, secondPoint) << endl;
00164   //   cout << "      min rms: " << computeMinRMS(firstPoint, secondPoint) << endl;
00165   if(firstPoint.deltaMean(secondPoint) <  computeAverageRMS(firstPoint, secondPoint) &&
00166      firstPoint.deltaRMS(secondPoint) < computeMinRMS(firstPoint, secondPoint)) {
00167     theValidity = true;
00168     //     cout << "    seed is valid!" << endl;
00169     return true;
00170   }  
00171   //   cout << "    seed is not valid!" << endl;
00172   theValidity = false;
00173   return false;
00174 }
00175 
00176 
00177 
00178   
00179 bool  DTOccupancyCluster::qualityCriterion(const DTOccupancyPoint& anotherPoint) {
00180   //   double minDistance = distance(anotherPoint);
00181 
00182   double minrms = 0;
00183   if(anotherPoint.rms() < averageRMS()) minrms = anotherPoint.rms();
00184   else minrms = averageRMS();
00185 
00186   //   cout << "  delta mean: " <<   fabs(averageMean() - anotherPoint.mean()) << endl;
00187   //   cout << "  average RMS: " << fabs(averageRMS() + anotherPoint.rms())/2. << endl;
00188   //   cout << "  delta rms: " <<  fabs(averageRMS() - anotherPoint.rms()) << endl;
00189   //   cout << " min rms: " << minrms << endl;
00190   //   if(fabs(averageMean() - anotherPoint.mean()) < fabs(averageRMS() + anotherPoint.rms())/2. &&
00191   //      fabs(averageRMS() - anotherPoint.rms()) < minrms) {
00192 
00193   if(fabs(averageMean() - anotherPoint.mean()) < averageRMS() &&
00194      fabs(averageRMS() - anotherPoint.rms()) < 2*minrms/3.) {
00195     theValidity = true;
00196     //     cout << "    point is valid!" << endl;
00197     return true;
00198   }  
00199   //   cout << "    point is not valid!" << endl;
00200   theValidity = false;
00201   return false;
00202 
00203   //   if(thePoints.size() == 1) return false;
00204   //   // Get the minimum distance from the point
00205   //   double minDistance = distance(anotherPoint);
00206   //   cout << "     min distance from cluster is: " << minDistance << endl;
00207   //   cout << "     radius is: " << radius << endl;
00208   //   // compare the minimum distance with the radius;
00209   //   if(minDistance < 2*radius) return true;
00210   //   else return false;
00211 }
00212 
00213 
00214 
00215 
00216 void DTOccupancyCluster::computeRadius() {
00217   double radius_squared = 0;
00218   for(vector<DTOccupancyPoint>::const_iterator pt_i = thePoints.begin();
00219       pt_i != thePoints.end(); ++pt_i) {
00220     for(vector<DTOccupancyPoint>::const_iterator pt_j = thePoints.begin();
00221         pt_j != thePoints.end(); ++pt_j) {
00222       radius_squared += TMath::Power(pt_i->distance(*pt_j),2);
00223     }
00224   }
00225   radius_squared = radius_squared/(2*TMath::Power(thePoints.size()+1,2));
00226   radius = sqrt(radius_squared);
00227   //   cout << "    the new cluster radius is: " << radius << endl;
00228 }
00229 
00230 
00231 
00232 int DTOccupancyCluster::nPoints() const {
00233   return thePoints.size();
00234 }
00235 
00236 
00237 set<DTLayerId> DTOccupancyCluster::getLayerIDs() const {
00238   set<DTLayerId> ret;
00239   for(vector<DTOccupancyPoint>::const_iterator point = thePoints.begin();
00240       point != thePoints.end(); ++point) {
00241     ret.insert((*point).layerId());
00242   }
00243   return ret;
00244 }
00245 
00246 
00247 bool clusterIsLessThan(const DTOccupancyCluster& clusterOne, const DTOccupancyCluster& clusterTwo) {
00248   if(clusterTwo.nPoints() == 1 && clusterOne.nPoints() != 1) {
00249     return true;
00250   }
00251   if(clusterTwo.nPoints() != 1 && clusterOne.nPoints() == 1) {
00252     return false;
00253   }
00254 
00255   if(clusterOne.nPoints() > clusterTwo.nPoints()) {
00256     return true;
00257   } else if(clusterOne.nPoints() < clusterTwo.nPoints()) {
00258     return false;
00259   } else {
00260     if(fabs(clusterOne.averageRMS() - sqrt(clusterOne.averageMean())) <
00261        fabs(clusterTwo.averageRMS() - sqrt(clusterTwo.averageMean()))) {
00262       return true;
00263     }
00264   }
00265   return false;
00266 
00267 }
00268 
00269