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