Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-05-01 22:51:49

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