File indexing completed on 2024-04-06 12:07:05
0001
0002
0003
0004
0005
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
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
0055 thePoints.push_back(singlePoint);
0056 }
0057
0058 DTOccupancyCluster::~DTOccupancyCluster() {}
0059
0060
0061 bool DTOccupancyCluster::isValid() const { return theValidity; }
0062
0063
0064
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
0072 computeRadius();
0073
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
0088
0089 double DTOccupancyCluster::distance(const DTOccupancyPoint& point) const {
0090 double dist = 99999999;
0091
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 }