Reference Manual

Go to the documentation of this file.
00001 
00002 #include "DQM/CSCMonitorModule/interface/CSCDQM_StripClusterFinder.h"
00003 
00004 namespace cscdqm {
00005 
00006 StripClusterFinder::StripClusterFinder(int l, int s, int cf, int st)
00007 {
00008 //
00009 // Options
00010 //
00011 //      fOpt = new CalibOptions();
00012         LayerNmb = l;
00013         TimeSliceNmb = s;
00014         StripNmb = cf*16;
00015 }
00016 void StripClusterFinder::DoAction(int LayerId, float *Cathodes)
00017 {
00018   int TimeId,StripId;
00019   this->LId=LayerId;
00020   MEStripClusters.clear();
00021   StripClusterFitData PulseHeightMapTMP;
00022   
00023   thePulseHeightMap.clear();
00024   
00025   // fill 
00026   //===========================================================================
00027 
00028   for(StripId=0;StripId<StripNmb;StripId++) { 
00029     for(TimeId=0;TimeId<TimeSliceNmb;TimeId++) { 
00030       PulseHeightMapTMP.height_[TimeId]=*(Cathodes+StripNmb*(TimeSliceNmb*LayerId+TimeId)+StripId);
00031     }
00032     PulseHeightMapTMP.bx_=0;
00033     PulseHeightMapTMP.channel_=StripId; // was StripId
00034     thePulseHeightMap.push_back(PulseHeightMapTMP);
00035   }
00036   sort(thePulseHeightMap.begin(),thePulseHeightMap.end(),Sort());
00037   //===========================================================================
00038   
00039   if(!thePulseHeightMap.size()) return;
00040   SearchMax();
00041   SearchBorders();
00042   Match();
00043   RefindMax();
00044   /*
00045   int val;
00046   for(i=0;i<MEStripClusters.size();i++){
00047     val=MEStripClusters[i].LFTBNDStrip;
00048     MEStripClusters[i].LFTBNDStrip=thePulseHeightMap[val].channel_;
00049     val=MEStripClusters[i].IRTBNDStrip;
00050     MEStripClusters[i].IRTBNDStrip=thePulseHeightMap[val].channel_;
00051     for(j=0;j<MEStripClusters[i].localMax.size();j++){
00052       val=MEStripClusters[i].localMax[j].Strip;
00053       MEStripClusters[i].localMax[j].Strip=thePulseHeightMap[val].channel_;
00054     }
00055   }
00056   */
00057 
00058   float sumstrip;
00059   float sumtime;
00060   float sumheight;
00061 
00062   for(i=0;i<MEStripClusters.size();i++){
00063     MEStripClusters[i].ClusterPulseMapHeight.clear();
00064     for(j=0;j<thePulseHeightMap.size();j++){
00065       if(thePulseHeightMap[j].channel_ >= MEStripClusters[i].LFTBNDStrip && 
00066          thePulseHeightMap[j].channel_ <= MEStripClusters[i].IRTBNDStrip) 
00067         MEStripClusters[i].ClusterPulseMapHeight.push_back(thePulseHeightMap[j]);
00068     }
00069     sumstrip=0;
00070     sumtime=0;
00071     sumheight=0;
00072     for(uint32_t k=0;k<MEStripClusters[i].ClusterPulseMapHeight.size();k++){
00073       for(int l=0;l<16;l++){
00074         sumstrip+=
00075           MEStripClusters[i].ClusterPulseMapHeight[k].height_[l]*
00076           (MEStripClusters[i].ClusterPulseMapHeight[k].channel_+1);
00077         sumtime+=
00078           MEStripClusters[i].ClusterPulseMapHeight[k].height_[l]*(l+1);
00079         sumheight+=
00080           MEStripClusters[i].ClusterPulseMapHeight[k].height_[l];
00081       }
00082     }
00083     if (sumheight) {
00084         MEStripClusters[i].Mean[0]=sumstrip/sumheight;
00085         MEStripClusters[i].Mean[1]=sumtime/sumheight;
00086     }
00087   }
00088 //  printClusters();
00089   return;
00090 }
00091 
00092 void StripClusterFinder::SearchMax(void)
00093 {
00094   StripCluster tmpCluster;
00095   for(i=1;i<(thePulseHeightMap.size()-1);i++){
00096     if(thePulseHeightMap[i].channel_==63 || thePulseHeightMap[i].channel_==64) continue; 
00097     for(j=1;j<15;j++){
00098       if(thePulseHeightMap[i].height_[j]>thePulseHeightMap[i-1].height_[j  ] &&
00099          thePulseHeightMap[i].height_[j]>thePulseHeightMap[i+1].height_[j  ] &&
00100          thePulseHeightMap[i].height_[j]>thePulseHeightMap[i  ].height_[j-1] &&
00101          thePulseHeightMap[i].height_[j]>thePulseHeightMap[i  ].height_[j+1] &&
00102          thePulseHeightMap[i].height_[j]>thePulseHeightMap[i-1].height_[j-1] &&
00103          thePulseHeightMap[i].height_[j]>thePulseHeightMap[i-1].height_[j+1] &&
00104          thePulseHeightMap[i].height_[j]>thePulseHeightMap[i+1].height_[j-1] &&
00105          thePulseHeightMap[i].height_[j]>thePulseHeightMap[i+1].height_[j+1]) {
00106         
00107         tmpCluster.localMax.clear();
00108         localMaxTMP.Strip = i;
00109         localMaxTMP.Time = j;
00110         tmpCluster.localMax.push_back(localMaxTMP);
00111         tmpCluster.LayerId   = this->LId;
00112         tmpCluster.LFTBNDTime   = -100;
00113         tmpCluster.LFTBNDStrip  = -100;
00114         tmpCluster.IRTBNDTime   = -100;
00115         tmpCluster.IRTBNDStrip  = -100;      
00116         MEStripClusters.push_back(tmpCluster);
00117       }
00118     }
00119   }
00120   return;
00121 }
00122 void StripClusterFinder::SearchBorders(void)
00123 {
00124   uint32_t iS,iT,iL,jL,iR,jR;
00125   
00126   //              SEARCHING PARAMETERS OF THE CLASTERS (LEFT DOWN & RIGHT UP)
00127   
00128   for(i=0;i<MEStripClusters.size();i++){
00129     if(!MEStripClusters[i].localMax.size()) {
00130       std::cout<<"!!!Warning Cluster has'nt local Maxima"<<std::endl;
00131       continue;
00132     }
00133     iS=MEStripClusters[i].localMax[0].Strip;
00134     iT=MEStripClusters[i].localMax[0].Time;
00135     //              LEFT DOWN 
00136     // strip
00137     MEStripClusters[i].LFTBNDStrip=0;
00138     for(iL=iS-1;iL>0;iL--) {
00139       if(thePulseHeightMap[iL].channel_==64) {
00140         MEStripClusters[i].LFTBNDStrip=iL;
00141         break;
00142       }
00143       if(thePulseHeightMap[iL].height_[iT]==0.) {
00144         MEStripClusters[i].LFTBNDStrip=iL+1;
00145         break;
00146       }
00147     }
00148     //time
00149     MEStripClusters[i].LFTBNDTime=0;
00150     for(jL=iT-1;jL>0;jL--){         
00151       if(thePulseHeightMap[iS].height_[jL]==0.){
00152         MEStripClusters[i].LFTBNDTime=jL+1;
00153         break;
00154       }
00155     }
00156     //              RIGHT UP
00157     //strip
00158     MEStripClusters[i].IRTBNDStrip=thePulseHeightMap.size()-1;
00159     for(iR=iS+1;iR<thePulseHeightMap.size();iR++){
00160       if(thePulseHeightMap[iR].channel_==63){
00161         MEStripClusters[i].IRTBNDStrip=iR;
00162         break;
00163       }
00164       if(thePulseHeightMap[iR].height_[iT]==0.){
00165         MEStripClusters[i].IRTBNDStrip=iR-1;
00166         break;
00167       }
00168     }
00169     //time
00170     MEStripClusters[i].IRTBNDTime=15; 
00171     for(jR=iT+1;jR<16;jR++){
00172       if(thePulseHeightMap[iS].height_[jR]==0.){
00173         MEStripClusters[i].IRTBNDTime=jR-1;
00174         break;
00175       }
00176     }
00177   }
00178   return;
00179 }
00180     
00181 void StripClusterFinder::Match(void)
00182 {
00183   //              MATCHING THE OVERLAPING CLASTERS
00184   bool find2match;
00185   find2match=true;
00186   do {
00187     find2match=FindAndMatch();
00188   }while(find2match);
00189   
00190   return;
00191 }
00192 
00193 bool StripClusterFinder::FindAndMatch(void)
00194 {
00195   // Find clusters to match
00196   icstart=0; 
00197   for(ic1=icstart;ic1<MEStripClusters.size();ic1++){
00198     IC1MIN=MEStripClusters[ic1].LFTBNDStrip;
00199     IC1MAX=MEStripClusters[ic1].IRTBNDStrip;
00200     JC1MIN=MEStripClusters[ic1].LFTBNDTime;
00201     JC1MAX=MEStripClusters[ic1].IRTBNDTime;
00202     for(ic2=ic1+1;ic2<MEStripClusters.size();ic2++){
00203       IC2MIN=MEStripClusters[ic2].LFTBNDStrip;
00204       IC2MAX=MEStripClusters[ic2].IRTBNDStrip;
00205       JC2MIN=MEStripClusters[ic2].LFTBNDTime;
00206       JC2MAX=MEStripClusters[ic2].IRTBNDTime;
00207       if((IC2MIN>=IC1MIN && IC2MIN<=IC1MAX && 
00208           JC2MIN>=JC1MIN && JC2MIN<=JC1MAX    ) || 
00209          (IC2MIN>=IC1MIN && IC2MIN<=IC1MAX && 
00210           JC2MAX>=JC1MIN && JC2MAX<=JC1MAX    ) || 
00211          (IC2MAX>=IC1MIN && IC2MAX<=IC1MAX && 
00212           JC2MIN>=JC1MIN && JC2MIN<=JC1MAX    ) || 
00213          (IC2MAX>=IC1MIN && IC2MAX<=IC1MAX && 
00214           JC2MAX>=JC1MIN && JC2MAX<=JC1MAX    )    ){
00215         KillCluster();
00216         return true;
00217       }
00218       else {
00219         if((IC1MIN>=IC2MIN && IC1MIN<=IC2MAX && 
00220             JC1MIN>=JC2MIN && JC1MIN<=JC2MAX ) || 
00221            (IC1MIN>=IC2MIN && IC1MIN<=IC2MAX && 
00222             JC1MAX>=JC2MIN && JC1MAX<=JC2MAX ) || 
00223            (IC1MAX>=IC2MIN && IC1MAX<=IC2MAX && 
00224             JC1MIN>=JC2MIN && JC1MIN<=JC2MAX ) || 
00225            (IC1MAX>=IC2MIN && IC1MAX<=IC2MAX && 
00226             JC1MAX>=JC2MIN && JC1MAX<=JC2MAX )    ){
00227           KillCluster();
00228           return true;
00229         }
00230       }
00231     }
00232   }
00233   return false;
00234 }
00235 void StripClusterFinder::KillCluster(void)
00236 {
00237   // Match Clusters and kill one of clusters.
00238   if(IC1MIN<IC2MIN)
00239     MEStripClusters[ic1].LFTBNDStrip=IC1MIN;
00240   else
00241     MEStripClusters[ic1].LFTBNDStrip=IC2MIN;
00242   if(JC1MIN<JC2MIN)
00243     MEStripClusters[ic1].LFTBNDTime=JC1MIN;
00244   else
00245     MEStripClusters[ic1].LFTBNDTime=JC2MIN;
00246   if(IC1MAX>IC2MAX)
00247     MEStripClusters[ic1].IRTBNDStrip=IC1MAX;
00248   else
00249     MEStripClusters[ic1].IRTBNDStrip=IC2MAX;
00250   if(JC1MAX>JC2MAX)
00251     MEStripClusters[ic1].IRTBNDTime=JC1MAX;
00252   else
00253     MEStripClusters[ic1].IRTBNDTime=JC2MAX;
00254 
00255   MEStripClusters.erase(MEStripClusters.begin()+ic2);
00256   icstart=ic1;
00257 
00258   return;
00259 }
00260 void StripClusterFinder::RefindMax(void)
00261 {
00262   int iLS,iRS,iLT,iRT;
00263   int iS,jT;
00264   int ilocal;
00265   float GlobalMax;
00266   bool Erased;
00267   //             SEARCHING EXTREMUMS IN THE CLASTERS
00268   
00269   for(i=0;i<MEStripClusters.size();i++){
00270     MEStripClusters[i].localMax.clear();
00271     ilocal=0;
00272     iLS=MEStripClusters[i].LFTBNDStrip;
00273     iRS=MEStripClusters[i].IRTBNDStrip;
00274     iLT=MEStripClusters[i].LFTBNDTime;
00275     iRT=MEStripClusters[i].IRTBNDTime;
00276     
00277     /*
00278     for(iS=iLS+1;iS<=iRS-1;iS++){ 
00279       for(jT=iLT+1;jT<=iRT-1;jT++){
00280     */
00281     for(iS=iLS;iS<=iRS;iS++){ 
00282       if(thePulseHeightMap[iS].channel_==63 || thePulseHeightMap[iS].channel_==64) continue; 
00283       for(jT=iLT;jT<=iRT;jT++){
00284         if(iS==0 || jT ==0 || iS==79 || jT==7) continue;
00285         if(thePulseHeightMap[iS].height_[jT]>thePulseHeightMap[iS-1].height_[jT  ] &&
00286            thePulseHeightMap[iS].height_[jT]>thePulseHeightMap[iS+1].height_[jT  ] &&
00287            thePulseHeightMap[iS].height_[jT]>thePulseHeightMap[iS  ].height_[jT-1] &&
00288            thePulseHeightMap[iS].height_[jT]>thePulseHeightMap[iS  ].height_[jT+1] &&
00289            thePulseHeightMap[iS].height_[jT]>thePulseHeightMap[iS-1].height_[jT-1] &&
00290            thePulseHeightMap[iS].height_[jT]>thePulseHeightMap[iS-1].height_[jT+1] &&
00291            thePulseHeightMap[iS].height_[jT]>thePulseHeightMap[iS+1].height_[jT-1] &&
00292            thePulseHeightMap[iS].height_[jT]>thePulseHeightMap[iS+1].height_[jT+1]) {
00293           localMaxTMP.Strip=iS;
00294           localMaxTMP.Time=jT;
00295           MEStripClusters[i].localMax.push_back(localMaxTMP);
00296           ilocal++;
00297         }
00298       }
00299     }               
00300     // kill local maximums rellated to noise, maximums with pulse height less then 10% of Global max of clust.
00301     //fing Global Max
00302     GlobalMax=0;
00303     if(MEStripClusters[i].localMax.size()) {
00304       //std::cout << "Cluster: " << i << " Number of local maximums before erase: " 
00305       //                << MEStripClusters[i].localMax.size() << std::endl; 
00306       for(j=0;j<MEStripClusters[i].localMax.size();j++){
00307         iS=MEStripClusters[i].localMax[j].Strip;
00308         jT=MEStripClusters[i].localMax[j].Time;
00309         /*
00310           std::cout << "Current Max:" 
00311           << " " << iS
00312           << " " << jT
00313           << " " << thePulseHeightMap[iS].height_[jT] << std::endl;
00314         */
00315         if(thePulseHeightMap[iS].height_[jT]>GlobalMax)
00316           GlobalMax=thePulseHeightMap[iS].height_[jT];
00317       }
00318       GlobalMax=(float) (GlobalMax/10.); 
00319       //erase noise localMaximums
00320       do{
00321         Erased=false;
00322         for(j=0;j<MEStripClusters[i].localMax.size();j++){
00323           iS=MEStripClusters[i].localMax[j].Strip;
00324           jT=MEStripClusters[i].localMax[j].Time;
00325           if(thePulseHeightMap[iS].height_[jT]<GlobalMax){
00326             MEStripClusters[i].localMax.erase(MEStripClusters[i].localMax.begin()+j);
00327             Erased=true;
00328             break;
00329           }
00330         }
00331       } while(Erased);
00332 
00333       //debug outs
00334       //std::cout << "Cluster: " << i << " Number of local maximums: " 
00335       //        << MEStripClusters[i].localMax.size() << std::endl; 
00336       /*
00337       for(j=0;j<MEStripClusters[i].localMax.size();j++){
00338         iS=MEStripClusters[i].localMax[j].Strip;
00339         jT=MEStripClusters[i].localMax[j].Time;
00340         std::cout << "Local Max: " << j << " Strip: " << iS << " Time: " << jT 
00341                   << " Height: " << thePulseHeightMap[iS].height_[jT] 
00342                   << " Cut Value: " << GlobalMax << std::endl;
00343       }
00344       */
00345     }
00346   }
00347   return;
00348 }
00349 void StripClusterFinder::printClusters(void)
00350 {
00351   int iS,jT;
00352   std::cout << "====================================================================" << std::endl;     
00353   std::cout << "debug information from StripClusterFinder" << std::endl;        
00354   for(i=0;i<MEStripClusters.size();i++){
00355     if(!MEStripClusters[i].localMax.size()) continue;
00356     std::cout << " Cluster: " << i+1 
00357               << " Number of local Maximums " <<  MEStripClusters[i].localMax.size() << std::endl;
00358     for(j=0;j<MEStripClusters[i].localMax.size();j++){
00359       iS=MEStripClusters[i].localMax[j].Strip;
00360       jT=MEStripClusters[i].localMax[j].Time;
00361 
00362       //      std::cout << "Local Max: " << j << " Strip: " << iS << " Time: " << jT << std::endl;
00363       for(uint32_t k=0;k<MEStripClusters[i].ClusterPulseMapHeight.size();k++){
00364         if(MEStripClusters[i].ClusterPulseMapHeight[k].channel_==iS)
00365           std::cout << "Local Max: " << j+1 << " Strip: " << iS+1 << " Time: " << jT+1 
00366                     << " Height: " << MEStripClusters[i].ClusterPulseMapHeight[k].height_[jT] << std::endl;
00367       }
00368     }
00369     for(uint32_t k=0;k<MEStripClusters[i].ClusterPulseMapHeight.size();k++){
00370       std::cout << "Strip: " << MEStripClusters[i].ClusterPulseMapHeight[k].channel_+1;
00371       for(int l=0;l<16;l++)
00372         std::cout << " " << MEStripClusters[i].ClusterPulseMapHeight[k].height_[l];
00373       std::cout << std::endl;
00374     }
00375 
00376     std::cout << " Left  Top    corner strip: " << MEStripClusters[i].LFTBNDStrip+1 << " " 
00377               << " time: " << MEStripClusters[i].LFTBNDTime+1 << std::endl; 
00378     std::cout << " Right Bottom corner strip: " << MEStripClusters[i].IRTBNDStrip+1 << " " 
00379               << " time: " << MEStripClusters[i].IRTBNDTime+1 << std::endl; 
00380   }
00381  std::cout << "======================================================================" << std::endl;    
00382    return;
00383 }
00384 bool  StripClusterFinder::Sort::operator()(StripClusterFitData a , StripClusterFitData b) const
00385 {return  a.channel_ < b.channel_;}
00386 
00387 }
/data/refman/pasoursint/CMSSW_5_3_3/src/DQM/CSCMonitorModule/src/CSCDQM_StripClusterFinder.cc
