DTMonitorClient/src/DTOccupancyCluster.cc

0001
0002 /*
0003  *  See header file for a description of this class.
0004  *
0005  *  \author G. Cerminara - INFN Torino
0006  */
0007
0008 #include "DTOccupancyCluster.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010
0011 #include "TH2F.h"
0012 #include "TMath.h"
0013
0014 #include <iostream>
0015
0016 using namespace std;
0017 using namespace edm;
0018
0019 DTOccupancyCluster::DTOccupancyCluster(const DTOccupancyPoint& firstPoint, const DTOccupancyPoint& secondPoint)
0020     : radius(0), theMaxMean(-1.), theMaxRMS(-1.), theMeanSum(0.), theRMSSum(0.) {
0021   if (!qualityCriterion(firstPoint, secondPoint)) {
0022     theValidity = false;
0023   } else {
0024     // compute the cluster quantities
0025     thePoints.push_back(firstPoint);
0026     thePoints.push_back(secondPoint);
0027     if (firstPoint.mean() > secondPoint.mean())
0028       theMaxMean = firstPoint.mean();
0029     else
0030       theMaxMean = secondPoint.mean();
0031
0032     if (firstPoint.rms() > secondPoint.rms())
0033       theMaxRMS = firstPoint.rms();
0034     else
0035       theMaxRMS = secondPoint.rms();
0036     theMeanSum += firstPoint.mean();
0037     theRMSSum += firstPoint.rms();
0038
0039     theMeanSum += secondPoint.mean();
0040     theRMSSum += secondPoint.rms();
0041
0042     computeRadius();
0043   }
0044 }
0045
0046 DTOccupancyCluster::DTOccupancyCluster(const DTOccupancyPoint& singlePoint)
0047     : radius(0),
0048       theMaxMean(singlePoint.mean()),
0049       theMaxRMS(singlePoint.rms()),
0050       theMeanSum(singlePoint.mean()),
0051       theRMSSum(singlePoint.rms()) {
0052   theValidity = true;
0053
0054   // compute the cluster quantities
0055   thePoints.push_back(singlePoint);
0056 }
0057
0058 DTOccupancyCluster::~DTOccupancyCluster() {}
0059
0060 // Check if the cluster candidate satisfies the quality requirements
0061 bool DTOccupancyCluster::isValid() const { return theValidity; }
0062
0063 // Add a point to the cluster: returns false if the point does not satisfy the
0064 // quality requirement
0065 bool DTOccupancyCluster::addPoint(const DTOccupancyPoint& anotherPoint) {
0066   LogTrace("DTDQM|DTMonitorClient|DTOccupancyTest|DTOccupancyCluster")
0067       << "   Add a point to the cluster: mean: " << anotherPoint.mean() << " rms: " << anotherPoint.rms() << endl;
0068   if (qualityCriterion(anotherPoint)) {
0069     LogTrace("DTDQM|DTMonitorClient|DTOccupancyTest|DTOccupancyCluster") << "   point is valid" << endl;
0070     thePoints.push_back(anotherPoint);
0071     // Compute the new cluster size
0072     computeRadius();
0073     // compute the max mean and RMS
0074     if (anotherPoint.mean() > theMaxMean) {
0075       theMaxMean = anotherPoint.mean();
0076     }
0077     if (anotherPoint.rms() > theMaxRMS) {
0078       theMaxRMS = anotherPoint.rms();
0079     }
0080     theMeanSum += anotherPoint.mean();
0081     theRMSSum += anotherPoint.rms();
0082     return true;
0083   }
0084   return false;
0085 }
0086
0087 // Compute the distance of a single point from the cluster
0088 // (minimum distance with respect to the cluster points)
0089 double DTOccupancyCluster::distance(const DTOccupancyPoint& point) const {
0090   double dist = 99999999;
0091   // compute the minimum distance from a point
0092   for (vector<DTOccupancyPoint>::const_iterator pt = thePoints.begin(); pt != thePoints.end(); ++pt) {
0093     double distance = point.distance(*pt);
0094     if (distance < dist) {
0095       dist = distance;
0096     }
0097   }
0098   return dist;
0099 }
0100
0101 double DTOccupancyCluster::averageMean() const { return theMeanSum / (double)thePoints.size(); }
0102
0103 double DTOccupancyCluster::averageRMS() const { return theRMSSum / (double)thePoints.size(); }
0104
0105 double DTOccupancyCluster::maxMean() const { return theMaxMean; }
0106
0107 double DTOccupancyCluster::maxRMS() const { return theMaxRMS; }
0108
0109 TH2F* DTOccupancyCluster::getHisto(std::string histoName,
0110                                    int nBinsX,
0111                                    double minX,
0112                                    double maxX,
0113                                    int nBinsY,
0114                                    double minY,
0115                                    double maxY,
0116                                    int fillColor) const {
0117   TH2F* histo = new TH2F(histoName.c_str(), histoName.c_str(), nBinsX, minX, maxX, nBinsY, minY, maxY);
0118   histo->SetFillColor(fillColor);
0119   for (vector<DTOccupancyPoint>::const_iterator pt = thePoints.begin(); pt != thePoints.end(); ++pt) {
0120     histo->Fill((*pt).mean(), (*pt).rms());
0121   }
0122   return histo;
0123 }
0124
0125 bool DTOccupancyCluster::qualityCriterion(const DTOccupancyPoint& firstPoint, const DTOccupancyPoint& secondPoint) {
0126   if (firstPoint.deltaMean(secondPoint) < computeAverageRMS(firstPoint, secondPoint) &&
0127       firstPoint.deltaRMS(secondPoint) < computeMinRMS(firstPoint, secondPoint)) {
0128     theValidity = true;
0129
0130     return true;
0131   }
0132
0133   theValidity = false;
0134   return false;
0135 }
0136
0137 bool DTOccupancyCluster::qualityCriterion(const DTOccupancyPoint& anotherPoint) {
0138   double minrms = 0;
0139   if (anotherPoint.rms() < averageRMS())
0140     minrms = anotherPoint.rms();
0141   else
0142     minrms = averageRMS();
0143
0144   if (fabs(averageMean() - anotherPoint.mean()) < averageRMS() &&
0145       fabs(averageRMS() - anotherPoint.rms()) < 2 * minrms / 3.) {
0146     theValidity = true;
0147     return true;
0148   }
0149   theValidity = false;
0150   return false;
0151 }
0152
0153 void DTOccupancyCluster::computeRadius() {
0154   double radius_squared = 0;
0155   for (vector<DTOccupancyPoint>::const_iterator pt_i = thePoints.begin(); pt_i != thePoints.end(); ++pt_i) {
0156     for (vector<DTOccupancyPoint>::const_iterator pt_j = thePoints.begin(); pt_j != thePoints.end(); ++pt_j) {
0157       radius_squared += TMath::Power(pt_i->distance(*pt_j), 2);
0158     }
0159   }
0160   radius_squared = radius_squared / (2 * TMath::Power(thePoints.size() + 1, 2));
0161   radius = sqrt(radius_squared);
0162 }
0163
0164 int DTOccupancyCluster::nPoints() const { return thePoints.size(); }
0165
0166 set<DTLayerId> DTOccupancyCluster::getLayerIDs() const {
0167   set<DTLayerId> ret;
0168   for (vector<DTOccupancyPoint>::const_iterator point = thePoints.begin(); point != thePoints.end(); ++point) {
0169     ret.insert((*point).layerId());
0170   }
0171   return ret;
0172 }
0173
0174 bool clusterIsLessThan(const DTOccupancyCluster& clusterOne, const DTOccupancyCluster& clusterTwo) {
0175   if (clusterTwo.nPoints() == 1 && clusterOne.nPoints() != 1) {
0176     return true;
0177   }
0178   if (clusterTwo.nPoints() != 1 && clusterOne.nPoints() == 1) {
0179     return false;
0180   }
0181
0182   if (clusterOne.nPoints() > clusterTwo.nPoints()) {
0183     return true;
0184   } else if (clusterOne.nPoints() < clusterTwo.nPoints()) {
0185     return false;
0186   } else {
0187     if (fabs(clusterOne.averageRMS() - sqrt(clusterOne.averageMean())) <
0188         fabs(clusterTwo.averageRMS() - sqrt(clusterTwo.averageMean()))) {
0189       return true;
0190     }
0191   }
0192   return false;
0193 }