Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:06:57 

0001 
0002 #include "CSCDQM_StripClusterFinder.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 
0005 namespace cscdqm {
0006 
0007   StripClusterFinder::StripClusterFinder(int l, int s, int cf, int st, bool ME11) {
0008     //
0009     // Options
0010     //
0011     //  fOpt = new CalibOptions();
0012     LayerNmb = l;
0013     TimeSliceNmb = s;
0014     StripNmb = cf * 16;
0015     isME11 = ME11;
0016     if (cf == 7) {
0017       is7DCFEBs = true;
0018       isME11 = true;
0019     } else {
0020       is7DCFEBs = false;
0021     }
0022   }
0023   void StripClusterFinder::DoAction(int LayerId, float* Cathodes) {
0024     int TimeId, StripId;
0025     MEStripClusters.clear();
0026     StripClusterFitData PulseHeightMapTMP;
0027 
0028     thePulseHeightMap.clear();
0029 
0030     // fill
0031     //===========================================================================
0032 
0033     for (StripId = 0; StripId < StripNmb; StripId++) {
0034       for (TimeId = 0; TimeId < TimeSliceNmb; TimeId++) {
0035         PulseHeightMapTMP.height_[TimeId] = *(Cathodes + StripNmb * (TimeSliceNmb * LayerId + TimeId) + StripId);
0036       }
0037       PulseHeightMapTMP.bx_ = 0;
0038       PulseHeightMapTMP.channel_ = StripId;  // was StripId
0039       thePulseHeightMap.push_back(PulseHeightMapTMP);
0040     }
0041     sort(thePulseHeightMap.begin(), thePulseHeightMap.end(), Sort());
0042     //===========================================================================
0043 
0044     if (thePulseHeightMap.empty())
0045       return;
0046     SearchMax(LayerId);
0047     SearchBorders();
0048     Match();
0049     RefindMax();
0050     /*
0051   int val;
0052   for(i=0;i<MEStripClusters.size();i++){
0053     val=MEStripClusters[i].LFTBNDStrip;
0054     MEStripClusters[i].LFTBNDStrip=thePulseHeightMap[val].channel_;
0055     val=MEStripClusters[i].IRTBNDStrip;
0056     MEStripClusters[i].IRTBNDStrip=thePulseHeightMap[val].channel_;
0057     for(j=0;j<MEStripClusters[i].localMax.size();j++){
0058       val=MEStripClusters[i].localMax[j].Strip;
0059       MEStripClusters[i].localMax[j].Strip=thePulseHeightMap[val].channel_;
0060     }
0061   }
0062   */
0063 
0064     float sumstrip;
0065     float sumtime;
0066     float sumheight;
0067 
0068     for (uint32_t i = 0; i < MEStripClusters.size(); i++) {
0069       MEStripClusters[i].ClusterPulseMapHeight.clear();
0070       for (uint32_t j = 0; j < thePulseHeightMap.size(); j++) {
0071         if (thePulseHeightMap[j].channel_ >= MEStripClusters[i].LFTBNDStrip &&
0072             thePulseHeightMap[j].channel_ <= MEStripClusters[i].IRTBNDStrip)
0073           MEStripClusters[i].ClusterPulseMapHeight.push_back(thePulseHeightMap[j]);
0074       }
0075       sumstrip = 0;
0076       sumtime = 0;
0077       sumheight = 0;
0078       for (uint32_t k = 0; k < MEStripClusters[i].ClusterPulseMapHeight.size(); k++) {
0079         for (int l = 0; l < 16; l++) {
0080           sumstrip += MEStripClusters[i].ClusterPulseMapHeight[k].height_[l] *
0081                       (MEStripClusters[i].ClusterPulseMapHeight[k].channel_ + 1);
0082           sumtime += MEStripClusters[i].ClusterPulseMapHeight[k].height_[l] * (l + 1);
0083           sumheight += MEStripClusters[i].ClusterPulseMapHeight[k].height_[l];
0084         }
0085       }
0086       if (sumheight) {
0087         MEStripClusters[i].Mean[0] = sumstrip / sumheight;
0088         MEStripClusters[i].Mean[1] = sumtime / sumheight;
0089       }
0090     }
0091     //  printClusters();
0092     return;
0093   }
0094 
0095   void StripClusterFinder::SearchMax(int32_t layerId) {
0096     StripCluster tmpCluster;
0097     for (uint32_t i = 1; i < (thePulseHeightMap.size() - 1); i++) {
0098       if (isME11 && (thePulseHeightMap[i].channel_ == 63 || thePulseHeightMap[i].channel_ == 64))
0099         continue;
0100       for (uint32_t j = 1; j < 15; j++) {
0101         if (thePulseHeightMap[i].height_[j] > thePulseHeightMap[i - 1].height_[j] &&
0102             thePulseHeightMap[i].height_[j] > thePulseHeightMap[i + 1].height_[j] &&
0103             thePulseHeightMap[i].height_[j] > thePulseHeightMap[i].height_[j - 1] &&
0104             thePulseHeightMap[i].height_[j] > thePulseHeightMap[i].height_[j + 1] &&
0105             thePulseHeightMap[i].height_[j] > thePulseHeightMap[i - 1].height_[j - 1] &&
0106             thePulseHeightMap[i].height_[j] > thePulseHeightMap[i - 1].height_[j + 1] &&
0107             thePulseHeightMap[i].height_[j] > thePulseHeightMap[i + 1].height_[j - 1] &&
0108             thePulseHeightMap[i].height_[j] > thePulseHeightMap[i + 1].height_[j + 1]) {
0109           tmpCluster.localMax.clear();
0110           localMaxTMP.Strip = i;
0111           localMaxTMP.Time = j;
0112           tmpCluster.localMax.push_back(localMaxTMP);
0113           tmpCluster.LayerId = layerId;
0114           tmpCluster.LFTBNDTime = -100;
0115           tmpCluster.LFTBNDStrip = -100;
0116           tmpCluster.IRTBNDTime = -100;
0117           tmpCluster.IRTBNDStrip = -100;
0118           MEStripClusters.push_back(tmpCluster);
0119         }
0120       }
0121     }
0122     return;
0123   }
0124   void StripClusterFinder::SearchBorders(void) {
0125     uint32_t iS, iT, iL, jL, iR, jR;
0126 
0127     //              SEARCHING PARAMETERS OF THE CLASTERS (LEFT DOWN & RIGHT UP)
0128 
0129     for (uint32_t i = 0; i < MEStripClusters.size(); i++) {
0130       if (MEStripClusters[i].localMax.empty()) {
0131         edm::LogWarning("NoLocalMax") << "Cluster " << i << " has no local Maxima";
0132         continue;
0133       }
0134       iS = MEStripClusters[i].localMax[0].Strip;
0135       iT = MEStripClusters[i].localMax[0].Time;
0136       //              LEFT DOWN
0137       // strip
0138       MEStripClusters[i].LFTBNDStrip = 0;
0139       for (iL = iS - 1; iL > 0; iL--) {
0140         if (isME11 && (thePulseHeightMap[iL].channel_ == 64)) {
0141           MEStripClusters[i].LFTBNDStrip = iL;
0142           break;
0143         }
0144         if (thePulseHeightMap[iL].height_[iT] == 0.) {
0145           MEStripClusters[i].LFTBNDStrip = iL + 1;
0146           break;
0147         }
0148       }
0149       //time
0150       MEStripClusters[i].LFTBNDTime = 0;
0151       for (jL = iT - 1; jL > 0; jL--) {
0152         if (thePulseHeightMap[iS].height_[jL] == 0.) {
0153           MEStripClusters[i].LFTBNDTime = jL + 1;
0154           break;
0155         }
0156       }
0157       //              RIGHT UP
0158       //strip
0159       MEStripClusters[i].IRTBNDStrip = thePulseHeightMap.size() - 1;
0160       for (iR = iS + 1; iR < thePulseHeightMap.size(); iR++) {
0161         if (isME11 && (thePulseHeightMap[iR].channel_ == 63)) {
0162           MEStripClusters[i].IRTBNDStrip = iR;
0163           break;
0164         }
0165         if (thePulseHeightMap[iR].height_[iT] == 0.) {
0166           MEStripClusters[i].IRTBNDStrip = iR - 1;
0167           break;
0168         }
0169       }
0170       //time
0171       MEStripClusters[i].IRTBNDTime = 15;
0172       for (jR = iT + 1; jR < 16; jR++) {
0173         if (thePulseHeightMap[iS].height_[jR] == 0.) {
0174           MEStripClusters[i].IRTBNDTime = jR - 1;
0175           break;
0176         }
0177       }
0178     }
0179     return;
0180   }
0181 
0182   void StripClusterFinder::Match(void) {
0183     //              MATCHING THE OVERLAPING CLASTERS
0184     bool find2match = true;
0185     do {
0186       find2match = FindAndMatch();
0187     } while (find2match);
0188 
0189     return;
0190   }
0191 
0192   bool StripClusterFinder::FindAndMatch(void) {
0193     // Find clusters to match
0194     for (uint32_t ic1 = 0; ic1 < MEStripClusters.size(); ic1++) {
0195       C1 c1;
0196       c1.IC1MIN = MEStripClusters[ic1].LFTBNDStrip;
0197       c1.IC1MAX = MEStripClusters[ic1].IRTBNDStrip;
0198       c1.JC1MIN = MEStripClusters[ic1].LFTBNDTime;
0199       c1.JC1MAX = MEStripClusters[ic1].IRTBNDTime;
0200       for (uint32_t ic2 = ic1 + 1; ic2 < MEStripClusters.size(); ic2++) {
0201         C2 c2;
0202         c2.IC2MIN = MEStripClusters[ic2].LFTBNDStrip;
0203         c2.IC2MAX = MEStripClusters[ic2].IRTBNDStrip;
0204         c2.JC2MIN = MEStripClusters[ic2].LFTBNDTime;
0205         c2.JC2MAX = MEStripClusters[ic2].IRTBNDTime;
0206         if ((c2.IC2MIN >= c1.IC1MIN && c2.IC2MIN <= c1.IC1MAX && c2.JC2MIN >= c1.JC1MIN && c2.JC2MIN <= c1.JC1MAX) ||
0207             (c2.IC2MIN >= c1.IC1MIN && c2.IC2MIN <= c1.IC1MAX && c2.JC2MAX >= c1.JC1MIN && c2.JC2MAX <= c1.JC1MAX) ||
0208             (c2.IC2MAX >= c1.IC1MIN && c2.IC2MAX <= c1.IC1MAX && c2.JC2MIN >= c1.JC1MIN && c2.JC2MIN <= c1.JC1MAX) ||
0209             (c2.IC2MAX >= c1.IC1MIN && c2.IC2MAX <= c1.IC1MAX && c2.JC2MAX >= c1.JC1MIN && c2.JC2MAX <= c1.JC1MAX)) {
0210           KillCluster(ic1, ic2, c1, c2);
0211           return true;
0212         } else {
0213           if ((c1.IC1MIN >= c2.IC2MIN && c1.IC1MIN <= c2.IC2MAX && c1.JC1MIN >= c2.JC2MIN && c1.JC1MIN <= c2.JC2MAX) ||
0214               (c1.IC1MIN >= c2.IC2MIN && c1.IC1MIN <= c2.IC2MAX && c1.JC1MAX >= c2.JC2MIN && c1.JC1MAX <= c2.JC2MAX) ||
0215               (c1.IC1MAX >= c2.IC2MIN && c1.IC1MAX <= c2.IC2MAX && c1.JC1MIN >= c2.JC2MIN && c1.JC1MIN <= c2.JC2MAX) ||
0216               (c1.IC1MAX >= c2.IC2MIN && c1.IC1MAX <= c2.IC2MAX && c1.JC1MAX >= c2.JC2MIN && c1.JC1MAX <= c2.JC2MAX)) {
0217             KillCluster(ic1, ic2, c1, c2);
0218             return true;
0219           }
0220         }
0221       }
0222     }
0223     return false;
0224   }
0225   void StripClusterFinder::KillCluster(uint32_t ic1, uint32_t ic2, C1 const& c1, C2 const& c2) {
0226     // Match Clusters and kill one of clusters.
0227     if (c1.IC1MIN < c2.IC2MIN)
0228       MEStripClusters[ic1].LFTBNDStrip = c1.IC1MIN;
0229     else
0230       MEStripClusters[ic1].LFTBNDStrip = c2.IC2MIN;
0231     if (c1.JC1MIN < c2.JC2MIN)
0232       MEStripClusters[ic1].LFTBNDTime = c1.JC1MIN;
0233     else
0234       MEStripClusters[ic1].LFTBNDTime = c2.JC2MIN;
0235     if (c1.IC1MAX > c2.IC2MAX)
0236       MEStripClusters[ic1].IRTBNDStrip = c1.IC1MAX;
0237     else
0238       MEStripClusters[ic1].IRTBNDStrip = c2.IC2MAX;
0239     if (c1.JC1MAX > c2.JC2MAX)
0240       MEStripClusters[ic1].IRTBNDTime = c1.JC1MAX;
0241     else
0242       MEStripClusters[ic1].IRTBNDTime = c2.JC2MAX;
0243 
0244     MEStripClusters.erase(MEStripClusters.begin() + ic2);
0245     return;
0246   }
0247   void StripClusterFinder::RefindMax(void) {
0248     //             SEARCHING EXTREMUMS IN THE CLUSTERS
0249 
0250     for (uint32_t i = 0; i < MEStripClusters.size(); i++) {
0251       MEStripClusters[i].localMax.clear();
0252       int iLS = MEStripClusters[i].LFTBNDStrip;
0253       int iRS = MEStripClusters[i].IRTBNDStrip;
0254       int iLT = MEStripClusters[i].LFTBNDTime;
0255       int iRT = MEStripClusters[i].IRTBNDTime;
0256 
0257       for (int iS = iLS; iS <= iRS; iS++) {
0258         if (isME11 && (thePulseHeightMap[iS].channel_ == 63 || thePulseHeightMap[iS].channel_ == 64))
0259           continue;
0260         for (int jT = iLT; jT <= iRT; jT++) {
0261           if (iS == 0 || jT == 0 || (!is7DCFEBs && (iS == 79)) || (is7DCFEBs && (iS == 111)) || jT == 7)
0262             continue;
0263           if (thePulseHeightMap[iS].height_[jT] > thePulseHeightMap[iS - 1].height_[jT] &&
0264               thePulseHeightMap[iS].height_[jT] > thePulseHeightMap[iS + 1].height_[jT] &&
0265               thePulseHeightMap[iS].height_[jT] > thePulseHeightMap[iS].height_[jT - 1] &&
0266               thePulseHeightMap[iS].height_[jT] > thePulseHeightMap[iS].height_[jT + 1] &&
0267               thePulseHeightMap[iS].height_[jT] > thePulseHeightMap[iS - 1].height_[jT - 1] &&
0268               thePulseHeightMap[iS].height_[jT] > thePulseHeightMap[iS - 1].height_[jT + 1] &&
0269               thePulseHeightMap[iS].height_[jT] > thePulseHeightMap[iS + 1].height_[jT - 1] &&
0270               thePulseHeightMap[iS].height_[jT] > thePulseHeightMap[iS + 1].height_[jT + 1]) {
0271             localMaxTMP.Strip = iS;
0272             localMaxTMP.Time = jT;
0273             MEStripClusters[i].localMax.push_back(localMaxTMP);
0274           }
0275         }
0276       }
0277       // kill local maximums rellated to noise, maximums with pulse height less then 10% of Global max of clust.
0278       //fing Global Max
0279       float GlobalMax = 0;
0280       if (!MEStripClusters[i].localMax.empty()) {
0281         //std::cout << "Cluster: " << i << " Number of local maximums before erase: "
0282         //      << MEStripClusters[i].localMax.size() << std::endl;
0283         for (uint32_t j = 0; j < MEStripClusters[i].localMax.size(); j++) {
0284           int iS = MEStripClusters[i].localMax[j].Strip;
0285           int jT = MEStripClusters[i].localMax[j].Time;
0286           /*
0287       std::cout << "Current Max:" 
0288       << " " << iS
0289       << " " << jT
0290       << " " << thePulseHeightMap[iS].height_[jT] << std::endl;
0291     */
0292           if (thePulseHeightMap[iS].height_[jT] > GlobalMax)
0293             GlobalMax = thePulseHeightMap[iS].height_[jT];
0294         }
0295         GlobalMax = (float)(GlobalMax / 10.);
0296         //erase noise localMaximums
0297         bool Erased;
0298         do {
0299           Erased = false;
0300           for (uint32_t j = 0; j < MEStripClusters[i].localMax.size(); j++) {
0301             int iS = MEStripClusters[i].localMax[j].Strip;
0302             int jT = MEStripClusters[i].localMax[j].Time;
0303             if (thePulseHeightMap[iS].height_[jT] < GlobalMax) {
0304               MEStripClusters[i].localMax.erase(MEStripClusters[i].localMax.begin() + j);
0305               Erased = true;
0306               break;
0307             }
0308           }
0309         } while (Erased);
0310 
0311         //debug outs
0312         //std::cout << "Cluster: " << i << " Number of local maximums: "
0313         //  << MEStripClusters[i].localMax.size() << std::endl;
0314         /*
0315       for(j=0;j<MEStripClusters[i].localMax.size();j++){
0316     iS=MEStripClusters[i].localMax[j].Strip;
0317     jT=MEStripClusters[i].localMax[j].Time;
0318     std::cout << "Local Max: " << j << " Strip: " << iS << " Time: " << jT 
0319           << " Height: " << thePulseHeightMap[iS].height_[jT] 
0320           << " Cut Value: " << GlobalMax << std::endl;
0321       }
0322       */
0323       }
0324     }
0325     return;
0326   }
0327   void StripClusterFinder::printClusters(void) {
0328     int iS, jT;
0329     std::cout << "====================================================================" << std::endl;
0330     std::cout << "debug information from StripClusterFinder" << std::endl;
0331     for (uint32_t i = 0; i < MEStripClusters.size(); i++) {
0332       if (MEStripClusters[i].localMax.empty())
0333         continue;
0334       std::cout << " Cluster: " << i + 1 << " Number of local Maximums " << MEStripClusters[i].localMax.size()
0335                 << std::endl;
0336       for (uint32_t j = 0; j < MEStripClusters[i].localMax.size(); j++) {
0337         iS = MEStripClusters[i].localMax[j].Strip;
0338         jT = MEStripClusters[i].localMax[j].Time;
0339 
0340         //      std::cout << "Local Max: " << j << " Strip: " << iS << " Time: " << jT << std::endl;
0341         for (uint32_t k = 0; k < MEStripClusters[i].ClusterPulseMapHeight.size(); k++) {
0342           if (MEStripClusters[i].ClusterPulseMapHeight[k].channel_ == iS)
0343             std::cout << "Local Max: " << j + 1 << " Strip: " << iS + 1 << " Time: " << jT + 1
0344                       << " Height: " << MEStripClusters[i].ClusterPulseMapHeight[k].height_[jT] << std::endl;
0345         }
0346       }
0347       for (uint32_t k = 0; k < MEStripClusters[i].ClusterPulseMapHeight.size(); k++) {
0348         std::cout << "Strip: " << MEStripClusters[i].ClusterPulseMapHeight[k].channel_ + 1;
0349         for (int l = 0; l < 16; l++)
0350           std::cout << " " << MEStripClusters[i].ClusterPulseMapHeight[k].height_[l];
0351         std::cout << std::endl;
0352       }
0353 
0354       std::cout << " Left  Top    corner strip: " << MEStripClusters[i].LFTBNDStrip + 1 << " "
0355                 << " time: " << MEStripClusters[i].LFTBNDTime + 1 << std::endl;
0356       std::cout << " Right Bottom corner strip: " << MEStripClusters[i].IRTBNDStrip + 1 << " "
0357                 << " time: " << MEStripClusters[i].IRTBNDTime + 1 << std::endl;
0358     }
0359     std::cout << "======================================================================" << std::endl;
0360     return;
0361   }
0362   bool StripClusterFinder::Sort::operator()(const StripClusterFitData& a, const StripClusterFitData& b) const {
0363     return a.channel_ < b.channel_;
0364   }
0365 
0366 }  // namespace cscdqm
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!