StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEmcOldFinder.cxx
1 #include "StEmcOldFinder.h"
2 #include "StEvent.h"
3 #include "StEventTypes.h"
4 #include "TH2.h"
5 #include "TString.h"
6 #include "StEmcUtil/others/emcDetectorName.h"
7 #include "StEmcPreCluster.h"
8 //#define TTABLESORTER
9 #ifdef TTABLESORTER
10 #include <TTableSorter.h>
11 #else
12 #include "TArrayI.h"
13 #include "TMath.h"
14 #endif
15 #include "Stiostream.h"
16 #include "StMessMgr.h"
17 
18 ClassImp(StEmcOldFinder)
19 
20 //_____________________________________________________________________________
22 {
23  Float_t seed[] = {0.7, 0.1, 0.4, 0.4};
24  Float_t add
25  [] =
26  {
27  0.001, 0.001, 0.001, 0.001
28  };
29  Float_t all[] =
30  {
31  0.1,0.1,0.1,0.1
32  };
33  Int_t size[] =
34  {
35  4,1,5,5
36  };
37  for(Int_t i = 0; i < MAXDETBARREL; i++)
38  {
39  mEnergySeed[i] = seed[i];
40  mEnergyAdd[i] = add
41  [i];
42  mEnergyThresholdAll[i] = all[i];
43  mSizeMax[i] = size[i];
44  if(mEnergyThresholdAll[i]<mEnergySeed[i])
45  mEnergyThresholdAll[i]=mEnergySeed[i];
46  }
47 }
48 //_____________________________________________________________________________
49 StEmcOldFinder::~StEmcOldFinder()
50 {}
52 {
53  StEmcCollection *emc = event->emcCollection();
54  if(!emc)
55  return kFALSE;
56 
57  for(Int_t i = 0; i<MAXDETBARREL; i++)
58  {
59  Int_t det = i+1;
60  StDetectorId id = static_cast<StDetectorId>(det-1+kBarrelEmcTowerId);
61  StEmcDetector* detector=emc->detector(id);
62  if(detector)
63  findClustersInDetector(detector);
64  }
65  return kTRUE;
66 }
68 {
69  if(!detector)
70  return kFALSE;
71  Int_t det = (Int_t)(detector->detectorId()-kBarrelEmcTowerId)+1;
72  Int_t NM = detector->numberOfModules();
73  for(Int_t m = 1; m<=NM; m++)
74  findClustersInModule(det,detector->module(m));
75  LOG_DEBUG <<"Number of clusters found for detector "<<det<<" = "<<mColl[det-1]->GetSize()<<endm;
76 
77  return kTRUE;
78 }
80 {
81  if(!module)
82  return kFALSE;
83 
84  StSPtrVecEmcRawHit& hits=module->hits();
85 
86  mFirst = 0;
87  mLast = hits.size();
88  mNH = mLast - mFirst;
89  if(mNH<=0)
90  return kFALSE;
91 
92  mEnergy.Set(mNH); // set dimensions of the arrays for this module
93  mEW.Set(mNH);
94  mSW.Set(mNH);
95  mEnergy.Reset(); // reset arrays
96  mEW.Reset();
97  mSW.Reset();
98  mUsed.Set(mNH); // array with flags of mUsed hits in the module
99  mUsed.Reset();
100 
101  Int_t ih;
102  Int_t jh = mFirst;
103  for(ih=mFirst; ih<mLast; ih++) // fill arrays
104  {
105  jh=ih-mFirst;
106  Int_t cal = hits[ih]->calibrationType();
107  if(cal<128) // ONLY hits with calibrationType <128 should be mUsed for cluster finding...
108  {
109  mEnergy[jh] = hits[ih]->energy();
110  mEW[jh]=hits[ih]->eta();
111  mSW[jh]=abs(hits[ih]->sub());
112  }
113  else
114  {
115  mEnergy[jh] = 0;
116  mEW[jh]=hits[ih]->eta();
117  mSW[jh]=abs(hits[ih]->sub());
118  }
119  }
120 
121  // only 1 hit in this module
122  if(mNH == 1)
123  {
124  if(mEnergy[0] > mEnergyThresholdAll[det-1])
125  {
126  StEmcPreCluster *cl = mColl[det-1]->newCluster();
127  cl->addHit(hits[0]);
128  cl->update();
129  mNHit = 0;
130  }
131  return kTRUE;
132  }
133 #ifdef TTABLESORTER
134  TTableSorter index(mEnergy.GetArray(), mNH); // The mLast element is biggest
135 #else
136  TArrayI index(mNH);
137  TMath::Sort(mNH,mEnergy.GetArray(),index.GetArray(),0);
138 #endif
139  mHitsId.Set(10);
140  mHitsId.Reset();
141  Int_t i, ii,l;
142  Int_t loop = (det==BTOW)?2:1;
143  Float_t eClW=0;
144 
145  for(i=mNH-1; i>=0; i--) //loop from the more energetic hit to the less one
146  {
147 #ifdef TTABLESORTER
148  Int_t j = index.GetIndex(i); //get index of the hit
149 #else
150  Int_t j = index[i];
151 #endif
152  if(mEnergy[j] < mEnergySeed[det-1])
153  break; //if the hit is below threshold for cluster find -> break
154 
155  if(mUsed[j] == 0) // test if this hit is not mUsed
156  {
157  mHitsId[0]=j; // First hit in cluster
158  mNHit=1;
159  mUsed[j]=1;
160  mKeyEta=0;
161  mKeyPhi=0;
162  eClW = mEnergy[j]; // 19-apr
163  if(mSizeMax[det-1] == 1)
164  goto TESTClUSTER; // cluster = hit
165 
166  for(l=0; l<loop; l++)
167  {
168  for(ii=i-1; ii>=0; ii--)
169  {
170 #ifdef TTABLESORTER
171  int jj = index.GetIndex(ii);
172 #else
173  int jj = index[ii];
174 #endif
175  if(mEnergy[jj] < mEnergyAdd[det-1])
176  break; // 19-apr
177  if(mUsed[jj] == 0)
178  {
179  if(!testOnNeighbor(det, jj))
180  {
181  mUsed[jj]=1;
182  mHitsId[mNHit]=jj;
183  mNHit++;
184  eClW += mEnergy[jj]; // 19-apr
185  if(mNHit == mSizeMax[det-1])
186  goto TESTClUSTER;
187  }
188  }
189  }
190  }
191 TESTClUSTER:
192  if(mNHit>0)
193  {
194  if(eClW > mEnergyThresholdAll[det-1])
195  {
196  StEmcPreCluster *cl = mColl[det-1]->newCluster();
197  for(Int_t H = 0; H<mNHit; H++)
198  cl->addHit(hits[mHitsId[H]+mFirst]);
199  cl->update();
200  mNHit = 0;
201  }
202  else
203  {
204  Int_t jj;
205  for(Int_t i1=0; i1<mNHit; i1++)
206  {
207  jj = mHitsId[i1];
208  mUsed[jj] = 0;
209  }
210  }
211  }
212  }
213  }
214  return kTRUE;
215 }
216 Bool_t StEmcOldFinder::testOnNeighbor(Int_t det, Int_t jn)
217 {
218  mOverlapFlag = 0;
219 
220  // find SMD-eta clusters
221  if(det==BSMDE)
222  {
223  if(mNHit == 1)
224  {
225  mEtaFirst=mEW[mHitsId[0]];
226  mEtaLast=mEtaFirst;
227  mEnergyFirst=mEnergy[mHitsId[0]];
228  mEnergyLast=mEnergyFirst;
229  }
230 
231  if (mEtaFirst-mEW[jn] == 1)
232  {
233  if (mEnergy[jn]<mEnergyFirst)
234  {
235  mEtaFirst=mEW[jn];
236  mEnergyFirst=mEnergy[jn];
237  return kFALSE;
238  }
239  else
240  mOverlapFlag++;
241 
242  }
243  else if (mEW[jn]-mEtaLast == 1)
244  {
245  if (mEnergy[jn]<mEnergyLast)
246  {
247  mEtaLast=mEW[jn];
248  mEnergyLast=mEnergy[jn];
249  return kFALSE;
250  }
251  else
252  mOverlapFlag++;
253  }
254  else
255  return kTRUE;
256  }
257 
258  // find SMD-phi clusters
259  else if(det==BSMDP)
260  {
261  if(mNHit == 1)
262  {
263  mEtaSeed=mEW[mHitsId[0]];
264  mPhiFirst=mSW[mHitsId[0]];
265  mPhiLast=mPhiFirst;
266  mEnergyFirst=mEnergy[mHitsId[0]];
267  mEnergyLast=mEnergyFirst;
268  }
269 
270  if(mEtaSeed == mEW[jn]) // Same eta bin
271  {
272  if(mPhiFirst-mSW[jn] == 1)
273  {
274  if (mEnergy[jn]<mEnergyFirst)
275  {
276  mPhiFirst=mSW[jn];
277  mEnergyFirst=mEnergy[jn];
278  return kFALSE;
279  }
280  else
281  mOverlapFlag++;
282  }
283  else if(mSW[jn]-mPhiLast == 1)
284  {
285  if (mEnergy[jn]<mEnergyLast)
286  {
287  mPhiLast=mSW[jn];
288  mEnergyLast=mEnergy[jn];
289  return kFALSE;
290  }
291  else
292  mOverlapFlag++;
293  }
294  else
295  return kTRUE;
296  }
297  }
298 
299  // find tower clusters
300  else
301  {
302  // The vector mNb can be used for further studies/checks.
303  Int_t mNb[5] = {0,0,0,0,0};
304  Int_t js = mHitsId[0];
305 
306  if (mEW[js] == mEW[jn] && abs(mSW[js]-mSW[jn]) == 1 && mNb[2] == 0)
307  {
308  mNb[2] = 1;
309  return kFALSE;
310  }
311  else if(mSW[js] == mSW[jn] && abs(mEW[js]-mEW[jn]) == 1)
312  {
313  if ((mEW[js]-mEW[jn]) == -1 && mNb[0] == 0)
314  {
315  mNb[0] = 1;
316  return kFALSE;
317  }
318  else if ((mEW[js]-mEW[jn]) == 1 && mNb[1] == 0)
319  {
320  mNb[1] = 1;
321  return kFALSE;
322  }
323  else
324  return kTRUE;
325  }
326  else if(abs(mSW[js]-mSW[jn]) == 1 && abs(mEW[js]-mEW[jn]) == 1)
327  {
328  if ((mEW[js]-mEW[jn]) == -1 && mNb[3] == 0)
329  {
330  mNb[3] = 1;
331  return kFALSE;
332  }
333  else if ((mEW[js]-mEW[jn]) == 1 && mNb[4] == 0)
334  {
335  mNb[4] = 1;
336  return kFALSE;
337  }
338  else
339  return kTRUE;
340  }
341  }
342  return kTRUE;
343 }
Bool_t findClustersInModule(Int_t, StEmcModule *)
finds clusters in a BEMC module
StEmcPreCluster * newCluster()
creates a new cluster in the collection. Returns its pointer
void addHit(StEmcRawHit *)
add a hit to the cluster
Bool_t testOnNeighbor(Int_t, Int_t)
test for hits in neighbor strips
void update()
updates cluster information. Calculates eta,phi, ernergy, etc from the hits added to the cluster ...
virtual Bool_t findClusters(StEvent *)
finds clusters in a StEvent object
Bool_t findClustersInDetector(StEmcDetector *)
finds clusters in a given detector