StRoot: StPreEclMaker/StEmcOldFinder.cxx Source File

00001 #include "StEmcOldFinder.h"
00002 #include "StEvent.h"
00003 #include "StEventTypes.h"
00004 #include "TH2.h"
00005 #include "TString.h"
00006 #include "StEmcUtil/others/emcDetectorName.h"
00007 #include "StEmcPreCluster.h"
00008 #include <TTableSorter.h>
00009 #include "Stiostream.h"
00010 #include "StMessMgr.h"
00011 
00012 ClassImp(StEmcOldFinder)
00013 
00014 //_____________________________________________________________________________
00015 StEmcOldFinder::StEmcOldFinder():StEmcVirtualFinder()
00016 {
00017     Float_t seed[] = {0.7, 0.1, 0.4, 0.4};
00018     Float_t add
00019         []  =
00020             {
00021                 0.001, 0.001, 0.001, 0.001
00022             };
00023     Float_t all[]  =
00024         {
00025             0.1,0.1,0.1,0.1
00026         };
00027     Int_t   size[] =
00028         {
00029             4,1,5,5
00030         };
00031     for(Int_t i = 0; i < MAXDETBARREL; i++)
00032     {
00033         mEnergySeed[i] = seed[i];
00034         mEnergyAdd[i] = add
00035                             [i];
00036         mEnergyThresholdAll[i] = all[i];
00037         mSizeMax[i] = size[i];
00038         if(mEnergyThresholdAll[i]<mEnergySeed[i])
00039             mEnergyThresholdAll[i]=mEnergySeed[i];
00040     }
00041 }
00042 //_____________________________________________________________________________
00043 StEmcOldFinder::~StEmcOldFinder()
00044 {}
00045 Bool_t StEmcOldFinder::findClusters(StEvent* event)
00046 {
00047     StEmcCollection *emc = event->emcCollection();
00048     if(!emc)
00049         return kFALSE;
00050 
00051     for(Int_t i = 0; i<MAXDETBARREL; i++)
00052     {
00053         Int_t det = i+1;
00054         StDetectorId id = static_cast<StDetectorId>(det-1+kBarrelEmcTowerId);
00055         StEmcDetector* detector=emc->detector(id);
00056         if(detector)
00057             findClustersInDetector(detector);
00058     }
00059     return kTRUE;
00060 }
00061 Bool_t StEmcOldFinder::findClustersInDetector(StEmcDetector* detector)
00062 {
00063     if(!detector)
00064         return kFALSE;
00065     Int_t det = (Int_t)(detector->detectorId()-kBarrelEmcTowerId)+1;
00066     Int_t NM = detector->numberOfModules();
00067     for(Int_t m = 1; m<=NM; m++)
00068         findClustersInModule(det,detector->module(m));
00069     LOG_DEBUG <<"Number of clusters found for detector "<<det<<" = "<<mColl[det-1]->GetSize()<<endm;
00070 
00071     return kTRUE;
00072 }
00073 Bool_t StEmcOldFinder::findClustersInModule(Int_t det, StEmcModule* module)
00074 {
00075     if(!module)
00076         return kFALSE;
00077 
00078     StSPtrVecEmcRawHit& hits=module->hits();
00079 
00080     mFirst = 0;
00081     mLast  = hits.size();
00082     mNH    = mLast - mFirst;
00083     if(mNH<=0)
00084         return kFALSE;
00085 
00086     mEnergy.Set(mNH);   // set dimensions of the arrays for this module
00087     mEW.Set(mNH);
00088     mSW.Set(mNH);
00089     mEnergy.Reset();   // reset arrays
00090     mEW.Reset();
00091     mSW.Reset();
00092     mUsed.Set(mNH);   // array with flags of mUsed hits in the module
00093     mUsed.Reset();
00094 
00095     Int_t  ih;
00096     Int_t  jh = mFirst;
00097     for(ih=mFirst; ih<mLast; ih++) // fill arrays
00098     {
00099         jh=ih-mFirst;
00100         Int_t cal = hits[ih]->calibrationType();
00101         if(cal<128) // ONLY hits with calibrationType <128 should be mUsed for cluster finding...
00102         {
00103             mEnergy[jh] = hits[ih]->energy();
00104             mEW[jh]=hits[ih]->eta();
00105             mSW[jh]=abs(hits[ih]->sub());
00106         }
00107         else
00108         {
00109             mEnergy[jh] = 0;
00110             mEW[jh]=hits[ih]->eta();
00111             mSW[jh]=abs(hits[ih]->sub());
00112         }
00113     }
00114 
00115     // only 1 hit in this module
00116     if(mNH == 1)
00117     {
00118         if(mEnergy[0] > mEnergyThresholdAll[det-1])
00119         {
00120             StEmcPreCluster *cl = mColl[det-1]->newCluster();
00121             cl->addHit(hits[0]);
00122             cl->update();
00123             mNHit = 0;
00124         }
00125         return kTRUE;
00126     }
00127 
00128     TTableSorter index(mEnergy.GetArray(), mNH); // The mLast element is biggest
00129     mHitsId.Set(10);
00130     mHitsId.Reset();
00131 
00132     Int_t i, ii,l;
00133     Int_t loop = (det==BTOW)?2:1;
00134     Float_t eClW=0;
00135 
00136     for(i=mNH-1; i>=0; i--) //loop from the more energetic hit to the less one
00137     {
00138         Int_t j = index.GetIndex(i); //get index of the hit
00139         if(mEnergy[j] < mEnergySeed[det-1])
00140             break; //if the hit is below threshold for cluster find -> break
00141 
00142         if(mUsed[j] == 0) // test if this hit is not mUsed
00143         {
00144             mHitsId[0]=j; // First hit in cluster
00145             mNHit=1;
00146             mUsed[j]=1;
00147             mKeyEta=0;
00148             mKeyPhi=0;
00149             eClW  = mEnergy[j]; // 19-apr
00150             if(mSizeMax[det-1] == 1)
00151                 goto TESTClUSTER; // cluster = hit
00152 
00153             for(l=0; l<loop; l++)
00154             {
00155                 for(ii=i-1; ii>=0; ii--)
00156                 {
00157                     int jj = index.GetIndex(ii);
00158                     if(mEnergy[jj] < mEnergyAdd[det-1])
00159                         break; // 19-apr
00160                     if(mUsed[jj] == 0)
00161                     {
00162                         if(!testOnNeighbor(det, jj))
00163                         {
00164                             mUsed[jj]=1;
00165                             mHitsId[mNHit]=jj;
00166                             mNHit++;
00167                             eClW += mEnergy[jj]; // 19-apr
00168                             if(mNHit == mSizeMax[det-1])
00169                                 goto TESTClUSTER;
00170                         }
00171                     }
00172                 }
00173             }
00174 TESTClUSTER:
00175             if(mNHit>0)
00176             {
00177                 if(eClW > mEnergyThresholdAll[det-1])
00178                 {
00179                     StEmcPreCluster *cl = mColl[det-1]->newCluster();
00180                     for(Int_t H = 0; H<mNHit; H++)
00181                         cl->addHit(hits[mHitsId[H]+mFirst]);
00182                     cl->update();
00183                     mNHit = 0;
00184                 }
00185                 else
00186                 {
00187                     Int_t jj;
00188                     for(Int_t i1=0; i1<mNHit; i1++)
00189                     {
00190                         jj = mHitsId[i1];
00191                         mUsed[jj] = 0;
00192                     }
00193                 }
00194             }
00195         }
00196     }
00197     return kTRUE;
00198 }
00199 Bool_t StEmcOldFinder::testOnNeighbor(Int_t det, Int_t jn)
00200 {
00201     mOverlapFlag = 0;
00202 
00203     // find SMD-eta clusters
00204     if(det==BSMDE)
00205     {
00206         if(mNHit == 1)
00207         {
00208             mEtaFirst=mEW[mHitsId[0]];
00209             mEtaLast=mEtaFirst;
00210             mEnergyFirst=mEnergy[mHitsId[0]];
00211             mEnergyLast=mEnergyFirst;
00212         }
00213 
00214         if (mEtaFirst-mEW[jn] == 1)
00215         {
00216             if (mEnergy[jn]<mEnergyFirst)
00217             {
00218                 mEtaFirst=mEW[jn];
00219                 mEnergyFirst=mEnergy[jn];
00220                 return kFALSE;
00221             }
00222             else
00223                 mOverlapFlag++;
00224 
00225         }
00226         else if (mEW[jn]-mEtaLast == 1)
00227         {
00228             if (mEnergy[jn]<mEnergyLast)
00229             {
00230                 mEtaLast=mEW[jn];
00231                 mEnergyLast=mEnergy[jn];
00232                 return kFALSE;
00233             }
00234             else
00235                 mOverlapFlag++;
00236         }
00237         else
00238             return kTRUE;
00239     }
00240 
00241     // find SMD-phi clusters
00242     else if(det==BSMDP)
00243     {
00244         if(mNHit == 1)
00245         {
00246             mEtaSeed=mEW[mHitsId[0]];
00247             mPhiFirst=mSW[mHitsId[0]];
00248             mPhiLast=mPhiFirst;
00249             mEnergyFirst=mEnergy[mHitsId[0]];
00250             mEnergyLast=mEnergyFirst;
00251         }
00252 
00253         if(mEtaSeed == mEW[jn])  // Same eta bin
00254         {
00255             if(mPhiFirst-mSW[jn] == 1)
00256             {
00257                 if (mEnergy[jn]<mEnergyFirst)
00258                 {
00259                     mPhiFirst=mSW[jn];
00260                     mEnergyFirst=mEnergy[jn];
00261                     return kFALSE;
00262                 }
00263                 else
00264                     mOverlapFlag++;
00265             }
00266             else if(mSW[jn]-mPhiLast == 1)
00267             {
00268                 if (mEnergy[jn]<mEnergyLast)
00269                 {
00270                     mPhiLast=mSW[jn];
00271                     mEnergyLast=mEnergy[jn];
00272                     return kFALSE;
00273                 }
00274                 else
00275                     mOverlapFlag++;
00276             }
00277             else
00278                 return kTRUE;
00279         }
00280     }
00281 
00282     // find tower clusters
00283     else
00284     {
00285         // The vector mNb can be used for further studies/checks.
00286         Int_t mNb[5] = {0,0,0,0,0};
00287         Int_t js = mHitsId[0];
00288 
00289         if     (mEW[js] == mEW[jn] && abs(mSW[js]-mSW[jn]) == 1 && mNb[2] == 0)
00290         {
00291             mNb[2] = 1;
00292             return kFALSE;
00293         }
00294         else if(mSW[js] == mSW[jn] && abs(mEW[js]-mEW[jn]) == 1)
00295         {
00296             if      ((mEW[js]-mEW[jn]) == -1 && mNb[0] == 0)
00297             {
00298                 mNb[0] = 1;
00299                 return kFALSE;
00300             }
00301             else if ((mEW[js]-mEW[jn]) == 1 && mNb[1] == 0)
00302             {
00303                 mNb[1] = 1;
00304                 return kFALSE;
00305             }
00306             else
00307                 return kTRUE;
00308         }
00309         else if(abs(mSW[js]-mSW[jn]) == 1 && abs(mEW[js]-mEW[jn]) == 1)
00310         {
00311             if      ((mEW[js]-mEW[jn]) == -1 && mNb[3] == 0)
00312             {
00313                 mNb[3] = 1;
00314                 return kFALSE;
00315             }
00316             else if ((mEW[js]-mEW[jn]) == 1 && mNb[4] == 0)
00317             {
00318                 mNb[4] = 1;
00319                 return kFALSE;
00320             }
00321             else
00322                 return kTRUE;
00323         }
00324     }
00325     return kTRUE;
00326 }