store also difference in local Y
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinder.cxx
CommitLineData
b0f5e3fc 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
04366a57 15////////////////////////////////////////////////////////////////////////////
16// //
17// Base Class used to find //
18// the reconstructed points for ITS //
19// See also AliITSClusterFinderSPD, AliITSClusterFinderSDD, //
20// AliITSClusterFinderSDD AliITSClusterFinderV2 //
21////////////////////////////////////////////////////////////////////////////
88cb7938 22
5d2c2f86 23#include "AliRun.h"
b0f5e3fc 24#include "AliITSClusterFinder.h"
04366a57 25#include "AliITSRecPoint.h"
7d62fb64 26#include "AliITSdigit.h"
27#include "AliITSDetTypeRec.h"
aacedc3e 28#include "AliITSMap.h"
1f3e997f 29#include "AliITSgeomTGeo.h"
5d2c2f86 30#include <TParticle.h>
6c8e94cf 31#include <TArrayI.h>
5d2c2f86 32#include "AliMC.h"
45a7397a 33#include "AliLog.h"
b0f5e3fc 34
fe7d86eb 35using std::endl;
36
b0f5e3fc 37ClassImp(AliITSClusterFinder)
38
5d2c2f86 39extern AliRun *gAlice;
40
9de0700b 41//----------------------------------------------------------------------
aacedc3e 42AliITSClusterFinder::AliITSClusterFinder():
43TObject(),
aacedc3e 44fModule(0),
45fDigits(0),
46fNdigits(0),
8ba39da9 47fDetTypeRec(0),
aacedc3e 48fClusters(0),
aacedc3e 49fMap(0),
5d2c2f86 50fNPeaks(-1),
51fNModules(AliITSgeomTGeo::GetNModules()),
d59a647b 52fEvent(0),
53fZmin(0),
54fZmax(0),
55fXmin(0),
6c8e94cf 56fXmax(0),
57fNClusters(0),
58fRawID2ClusID(0)
59{
9de0700b 60 // default cluster finder
aacedc3e 61 // Input:
62 // none.
63 // Output:
64 // none.
65 // Return:
66 // A default constructed AliITSCulsterFinder
40f5092a 67 for(Int_t i=0; i<2200; i++){
68 fNdet[i]=0;
69 fNlayer[i]=0;
70 }
9de0700b 71}
72//----------------------------------------------------------------------
8ba39da9 73AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp):
aacedc3e 74TObject(),
aacedc3e 75fModule(0),
76fDigits(0),
77fNdigits(0),
e56160b8 78fDetTypeRec(dettyp),
aacedc3e 79fClusters(0),
aacedc3e 80fMap(0),
5d2c2f86 81fNPeaks(-1),
82fNModules(AliITSgeomTGeo::GetNModules()),
d59a647b 83fEvent(0),
84fZmin(0),
85fZmax(0),
86fXmin(0),
6c8e94cf 87fXmax(0),
88fNClusters(0),
89fRawID2ClusID(0)
90{
5d2c2f86 91 // default cluster finder
aacedc3e 92 // Standard constructor for cluster finder
93 // Input:
94 // AliITSsegmentation *seg The segmentation class to be used
95 // AliITSresponse *res The response class to be used
96 // Output:
97 // none.
98 // Return:
99 // A Standard constructed AliITSCulsterFinder
40f5092a 100 for(Int_t i=0; i<2200; i++){
101 fNdet[i]=0;
102 fNlayer[i]=0;
103 }
aacedc3e 104}
105//----------------------------------------------------------------------
8ba39da9 106AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp,
107 TClonesArray *digits):
aacedc3e 108TObject(),
aacedc3e 109fModule(0),
110fDigits(digits),
111fNdigits(0),
8ba39da9 112fDetTypeRec(dettyp),
aacedc3e 113fClusters(0),
aacedc3e 114fMap(0),
5d2c2f86 115fNPeaks(-1),
116fNModules(AliITSgeomTGeo::GetNModules()),
d59a647b 117fEvent(0),
118fZmin(0),
119fZmax(0),
120fXmin(0),
6c8e94cf 121fXmax(0),
122fNClusters(0),
123fRawID2ClusID(0)
124{
5d2c2f86 125 // default cluster finder
aacedc3e 126 // Standard + cluster finder constructor
127 // Input:
128 // AliITSsegmentation *seg The segmentation class to be used
129 // AliITSresponse *res The response class to be used
130 // TClonesArray *digits Array of digits to be used
131 // Output:
132 // none.
133 // Return:
134 // A Standard constructed AliITSCulsterFinder
b0f5e3fc 135
40f5092a 136 fNdigits = fDigits->GetEntriesFast();
137 for(Int_t i=0; i<2200; i++){
138 fNdet[i]=0;
139 fNlayer[i]=0;
140 }
b0f5e3fc 141}
a971ed8e 142
04366a57 143//______________________________________________________________________
d59a647b 144AliITSClusterFinder::AliITSClusterFinder(const AliITSClusterFinder &source) :
145 TObject(source),
146 fModule(source.fModule),
147 fDigits(),
148 fNdigits(source.fNdigits),
149 fDetTypeRec(),
150 fClusters(),
151 fMap(),
152 fNPeaks(source.fNPeaks),
153 fNModules(source.fNModules),
154 fEvent(source.fEvent),
155 fZmin(source.fZmin),
156 fZmax(source.fZmax),
157 fXmin(source.fXmin),
6c8e94cf 158 fXmax(source.fXmax),
159 fNClusters(source.fNClusters),
160 fRawID2ClusID(source.fRawID2ClusID)
d59a647b 161{
04366a57 162 // Copy constructor
163 // Copies are not allowed. The method is protected to avoid misuse.
a971ed8e 164 AliError("Copy constructor not allowed\n");
04366a57 165}
a971ed8e 166
04366a57 167
168//______________________________________________________________________
e56160b8 169//AliITSClusterFinder& AliITSClusterFinder::operator=(const AliITSClusterFinder& /* source */){
04366a57 170 // Assignment operator
171 // Assignment is not allowed. The method is protected to avoid misuse.
e56160b8 172// Fatal("= operator","Assignment operator not allowed\n");
173// return *this;
174//}
04366a57 175
9de0700b 176//----------------------------------------------------------------------
177AliITSClusterFinder::~AliITSClusterFinder(){
178 // destructor cluster finder
aacedc3e 179 // Input:
180 // none.
181 // Output:
182 // none.
183 // Return:
184 // none.
b0f5e3fc 185
aacedc3e 186 if(fMap) {delete fMap;}
9de0700b 187 // Zero local pointers. Other classes own these pointers.
9de0700b 188 fMap = 0;
189 fDigits = 0;
190 fNdigits = 0;
9de0700b 191 fNPeaks = 0;
7d62fb64 192 fDetTypeRec = 0;
04366a57 193
7d62fb64 194}
b0f5e3fc 195//__________________________________________________________________________
7d62fb64 196void AliITSClusterFinder::InitGeometry(){
1f3e997f 197 //
198 // Initialisation of ITS geometry
199 //
200 Int_t mmax=AliITSgeomTGeo::GetNModules();
201 for (Int_t m=0; m<mmax; m++) {
202 Int_t lay,lad,det; AliITSgeomTGeo::GetModuleId(m,lay,lad,det);
1f3e997f 203 fNdet[m] = (lad-1)*AliITSgeomTGeo::GetNDetectors(lay) + (det-1);
04366a57 204 fNlayer[m] = lay-1;
205 }
b0f5e3fc 206}
04366a57 207
208
7d62fb64 209
04366a57 210
f8d9a5b8 211//______________________________________________________________________
aacedc3e 212Bool_t AliITSClusterFinder::IsNeighbor(TObjArray *digs,Int_t i,Int_t n[])const{
f8d9a5b8 213 // Locagical function which checks to see if digit i has a neighbor.
214 // If so, then it returns kTRUE and its neighbor index j.
215 // This routine checks if the digits are side by side or one before the
216 // other. Requires that the array of digits be in proper order.
217 // Returns kTRUE in the following cases.
218 // ji 0j if kdiagonal j0 0i
219 // 00 0i if kdiagonal 0i j0
220 // Inputs:
221 // TObjArray *digs Array to search for neighbors in
222 // Int_t i Index of digit for which we are searching for
223 // a neighbor of.
224 // Output:
225 // Int_t j[4] Index of one or more of the digits which is a
226 // neighbor of digit a index i.
227 // Return:
228 // Bool_t kTRUE if a neighbor was found kFALSE otherwise.
229 Int_t ix,jx,iz,jz,j;
230 const Bool_t kdiagonal=kFALSE;
231 Bool_t nei[4];
232
233 // No neighbors found if array empty.
234 if(digs->GetEntriesFast()<=0) return kFALSE;
235 // can not be a digit with first element or elements out or range
236 if(i<=0 || i>=digs->GetEntriesFast()) return kFALSE;
237
238 for(j=0;j<4;j++){n[j] = -1;nei[j]=kFALSE;}
239 ix = ((AliITSdigit*)(digs->At(i)))->GetCoord1();
240 iz = ((AliITSdigit*)(digs->At(i)))->GetCoord2();
241 for(j=0;j<i;j++){
aacedc3e 242 jx = ((AliITSdigit*)(digs->At(j)))->GetCoord1();
243 jz = ((AliITSdigit*)(digs->At(j)))->GetCoord2();
244 if(jx+1==ix && jz ==iz){n[0] = j;nei[0] = kTRUE;}
245 if(jx ==ix && jz+1==iz){n[1] = j;nei[1] = kTRUE;}
246 if(jx+1==ix && jz+1==iz){n[2] = j;nei[2] = kTRUE;}
247 if(jx+1==ix && jz-1==iz){n[3] = j;nei[3] = kTRUE;}
f8d9a5b8 248 } // end for k
249 if(nei[0]||nei[1]) return kTRUE;
250 if(kdiagonal&&(nei[2]||nei[3])) return kTRUE;
251 // no Neighbors found.
252 return kFALSE;
253}
aacedc3e 254
255//______________________________________________________________________
7d62fb64 256void AliITSClusterFinder::Print(ostream *os) const{
aacedc3e 257 //Standard output format for this class
258 // Inputs:
259 // ostream *os Output stream
260 // Output:
261 // ostream *os Output stream
262 // Return:
263 // none.
264
aacedc3e 265 *os << fModule<<",";
266 *os << fNdigits<<",";
aacedc3e 267 *os << fNPeaks<<endl;
268}
04366a57 269//______________________________________________________________________
7d62fb64 270void AliITSClusterFinder::Read(istream *is) {
aacedc3e 271 //Standard input for this class
272 // Inputs:
273 // istream *is Input stream
274 // Output:
275 // istream *is Input stream
276 // Return:
277 // none.
278
aacedc3e 279 *is >> fModule;
280 *is >> fNdigits;
aacedc3e 281 *is >> fNPeaks;
282}
283//______________________________________________________________________
284ostream &operator<<(ostream &os,AliITSClusterFinder &source){
285 // Standard output streaming function.
286 // Inputs:
287 // ostream *os Output stream
288 // AliITSClusterFinder &source Class to be printed
289 // Output:
290 // ostream *os Output stream
291 // Return:
292 // none.
293
294 source.Print(&os);
295 return os;
296}
297//______________________________________________________________________
298istream &operator>>(istream &is,AliITSClusterFinder &source){
299 // Standard output streaming function.
300 // Inputs:
301 // istream *is Input stream
302 // AliITSClusterFinder &source Class to be read in.
303 // Output:
304 // istream *is Input stream
305 // Return:
306 // none.
307
308 source.Read(&is);
309 return is;
310}
45a7397a 311
5d2c2f86 312//______________________________________________________________________
45a7397a 313void AliITSClusterFinder::CheckLabels2(Int_t lab[10])
314{
5d2c2f86 315 //------------------------------------------------------------
316 // Tries to find mother's labels
317 //------------------------------------------------------------
33c3c91a 318 AliRunLoader *rl = AliRunLoader::Instance();
c605e5f7 319 if(!rl) return;
5d2c2f86 320 TTree *trK=(TTree*)rl->TreeK();
45a7397a 321 if (!trK) return;
322 //
323 int labS[10];
324 Int_t nlabels = 0;
325 Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
326 for (Int_t i=0;i<10;i++) if (lab[i]>=0) labS[nlabels++] = lab[i];
327 if (nlabels==0) return;
328 //
329 float mom[10];
330 for (Int_t i=0;i<nlabels;i++) {
331 Int_t label = labS[i];
332 mom[i] = 0;
333 if (label>=ntracks) continue;
334 TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
335 mom[i] = part->P();
336 if (part->P() < 0.02) { // reduce soft particles from the same cluster
337 Int_t m=part->GetFirstMother();
338 if (m<0) continue; // primary
339 //
340 if (part->GetStatusCode()>0) continue;
341 //
342 // if the parent is within the same cluster, reassign the label to it
343 for (int j=0;j<nlabels;j++) if (labS[j]==m) { labS[i] = m; break; }
5d2c2f86 344 }
45a7397a 345 }
346 //
347 if (nlabels>3) { // only 3 labels are stored in cluster, sort in decreasing momentum
348 int ind[10],labSS[10];
349 TMath::Sort(nlabels,mom,ind);
350 for (int i=nlabels;i--;) labSS[i] = labS[i];
351 for (int i=0;i<nlabels;i++) labS[i] = labSS[ind[i]];
352 }
353 //
354 //compress labels -- if multi-times the same
355 for (Int_t i=0;i<10;i++) lab[i]=-2;
356 int nlabFin=0,j=0;
357 for (int i=0;i<nlabels;i++) {
358 for (j=0;j<nlabFin;j++) if (labS[i]==lab[j]) break; // the label already there
359 if (j==nlabFin) lab[nlabFin++] = labS[i];
5d2c2f86 360 }
45a7397a 361 //
5d2c2f86 362}
363
364//______________________________________________________________________
365void AliITSClusterFinder::AddLabel(Int_t lab[10], Int_t label) {
366 //add label to the cluster
33c3c91a 367 AliRunLoader *rl = AliRunLoader::Instance();
5d2c2f86 368 TTree *trK=(TTree*)rl->TreeK();
369 if(trK){
370 if(label<0) return; // In case of no label just exit
371
372 Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
373 if (label>ntracks) return;
374 for (Int_t i=0;i<10;i++){
375 // if (label<0) break;
376 if (lab[i]==label) break;
377 if (lab[i]<0) {
378 lab[i]= label;
379 break;
380 }
381 }
382 }
383}
384
385
386//______________________________________________________________________
387void AliITSClusterFinder::
388FindCluster(Int_t k,Int_t maxz,AliBin *bins,Int_t &n,Int_t *idx) {
389 //------------------------------------------------------------
390 // returns an array of indices of digits belonging to the cluster
391 // (needed when the segmentation is not regular)
392 //------------------------------------------------------------
393 if (n<200) idx[n++]=bins[k].GetIndex();
394 bins[k].Use();
395
396 if (bins[k-maxz].IsNotUsed()) FindCluster(k-maxz,maxz,bins,n,idx);
397 if (bins[k-1 ].IsNotUsed()) FindCluster(k-1 ,maxz,bins,n,idx);
398 if (bins[k+maxz].IsNotUsed()) FindCluster(k+maxz,maxz,bins,n,idx);
399 if (bins[k+1 ].IsNotUsed()) FindCluster(k+1 ,maxz,bins,n,idx);
400 /*
401 if (bins[k-maxz-1].IsNotUsed()) FindCluster(k-maxz-1,maxz,bins,n,idx);
402 if (bins[k-maxz+1].IsNotUsed()) FindCluster(k-maxz+1,maxz,bins,n,idx);
403 if (bins[k+maxz-1].IsNotUsed()) FindCluster(k+maxz-1,maxz,bins,n,idx);
404 if (bins[k+maxz+1].IsNotUsed()) FindCluster(k+maxz+1,maxz,bins,n,idx);
405 */
406}
407
408//______________________________________________________________________
409Bool_t AliITSClusterFinder::IsMaximum(Int_t k,Int_t max,const AliBin *bins) {
410 //------------------------------------------------------------
411 //is this a local maximum ?
412 //------------------------------------------------------------
413 UShort_t q=bins[k].GetQ();
414 if (q==1023) return kFALSE;
415 if (bins[k-max].GetQ() > q) return kFALSE;
416 if (bins[k-1 ].GetQ() > q) return kFALSE;
417 if (bins[k+max].GetQ() > q) return kFALSE;
418 if (bins[k+1 ].GetQ() > q) return kFALSE;
419 if (bins[k-max-1].GetQ() > q) return kFALSE;
420 if (bins[k+max-1].GetQ() > q) return kFALSE;
421 if (bins[k+max+1].GetQ() > q) return kFALSE;
422 if (bins[k-max+1].GetQ() > q) return kFALSE;
423 return kTRUE;
424}
425
426//______________________________________________________________________
427void AliITSClusterFinder::
428FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
429 //------------------------------------------------------------
430 //find local maxima
431 //------------------------------------------------------------
432 if (n<31)
433 if (IsMaximum(k,max,b)) {
434 idx[n]=k; msk[n]=(2<<n);
435 n++;
436 }
437 b[k].SetMask(0);
438 if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
439 if (b[k-1 ].GetMask()&1) FindPeaks(k-1 ,max,b,idx,msk,n);
440 if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
441 if (b[k+1 ].GetMask()&1) FindPeaks(k+1 ,max,b,idx,msk,n);
442}
443
444//______________________________________________________________________
445void AliITSClusterFinder::
446MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
447 //------------------------------------------------------------
448 //mark this peak
449 //------------------------------------------------------------
450 UShort_t q=bins[k].GetQ();
451
452 bins[k].SetMask(bins[k].GetMask()|m);
453
454 if (bins[k-max].GetQ() <= q)
455 if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
456 if (bins[k-1 ].GetQ() <= q)
457 if ((bins[k-1 ].GetMask()&m) == 0) MarkPeak(k-1 ,max,bins,m);
458 if (bins[k+max].GetQ() <= q)
459 if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
460 if (bins[k+1 ].GetQ() <= q)
461 if ((bins[k+1 ].GetMask()&m) == 0) MarkPeak(k+1 ,max,bins,m);
462}
463
464//______________________________________________________________________
465void AliITSClusterFinder::
466MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSRecPoint &c) {
467 //------------------------------------------------------------
468 //make cluster using digits of this peak
469 //------------------------------------------------------------
470 Float_t q=(Float_t)bins[k].GetQ();
471 Int_t i=k/max, j=k-i*max;
d59a647b 472 if(c.GetQ()<0.01){ // first entry in cluster
473 fXmin=i;
474 fXmax=i;
475 fZmin=j;
476 fZmax=j;
477 }else{ // check cluster extension
478 if(i<fXmin) fXmin=i;
479 if(i>fXmax) fXmax=i;
480 if(j<fZmin) fZmin=j;
481 if(j>fZmax) fZmax=j;
482 }
5d2c2f86 483 c.SetQ(c.GetQ()+q);
484 c.SetY(c.GetY()+i*q);
485 c.SetZ(c.GetZ()+j*q);
486 c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
487 c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);
488
489 bins[k].SetMask(0xFFFFFFFE);
6c8e94cf 490 if (fRawID2ClusID) { // RS: Register cluster id in raw words list
491 int rwid = bins[k].GetRawID();
492 if (fRawID2ClusID->GetSize()<=rwid) fRawID2ClusID->Set( (rwid+10)<<1 );
493 (*fRawID2ClusID)[rwid] = fNClusters+1; // RS: store clID+1 as a reference to the cluster
494 }
5d2c2f86 495 if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
496 if (bins[k-1 ].GetMask() == m) MakeCluster(k-1 ,max,bins,m,c);
497 if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
498 if (bins[k+1 ].GetMask() == m) MakeCluster(k+1 ,max,bins,m,c);
499}