don't sort clusters after local reco, do this in AliITSUTrackerGlo
[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSUClusterizer.cxx
CommitLineData
5e375bb4 1#include <TTree.h>
6cae87c5 2#include <TFile.h>
5e375bb4 3#include <TObjArray.h>
4bac12be 4#include <TParticle.h>
5e375bb4 5#include <TMath.h>
4bac12be 6#include "AliRun.h"
6f54924d 7#include "AliLog.h"
4bac12be 8#include "AliMC.h"
5e375bb4 9#include <AliITSUSegmentationPix.h>
b69620f8 10#include "AliITSUClusterizer.h"
5e375bb4 11#include "AliITSUGeomTGeo.h"
12#include "AliITSUSegmentationPix.h"
c7811d12 13#include "AliITSUDigitPix.h"
5e375bb4 14#include "AliITSURecoParam.h"
889b1493 15#include "AliITSUAux.h"
5e375bb4 16using namespace TMath;
889b1493 17using namespace AliITSUAux;
b69620f8 18
19ClassImp(AliITSUClusterizer)
20
21//______________________________________________________________________________
22AliITSUClusterizer::AliITSUClusterizer(Int_t initNRow)
23: fVolID(-1)
ee58ce21 24 ,fAllowDiagonalClusterization(kFALSE)
b69620f8 25 ,fSegm(0)
5e375bb4 26 ,fRecoParam(0)
b69620f8 27 ,fInputDigits(0)
28 ,fInputDigitsReadIndex(0)
889b1493 29 ,fLayerID(0)
4bac12be 30 ,fNLabels(0)
31 ,fRawData(kFALSE)
889b1493 32 ,fLorAngCorrection(0)
b69620f8 33 ,fOutputClusters(0)
34 ,fDigitFreelist(0)
35 ,fPartFreelist(0)
36 ,fCandFreelist(0)
37 ,fDigitFreelistBptrFirst(0)
38 ,fDigitFreelistBptrLast(0)
39 ,fPartFreelistBptr(0)
40 ,fCandFreelistBptr(0)
6cae87c5 41#ifdef _ClusterTopology_
6cae87c5 42 ,fTopology(0)
43 ,fMinCol(0)
44 ,fMinRow(0)
45#endif //_ClusterTopology_
b69620f8 46{
47 SetUniqueID(0);
48 // c-tor
49 SetNRow(initNRow);
6f54924d 50//
51#ifdef _ClusterTopology_
52 AliInfo("*********************************************************************");
53 AliInfo("ATTENTION: YOU ARE RUNNING IN SPECIAL MODE OF STORING CLUSTER PATTERN");
54 AliInfo("*********************************************************************");
55#endif //_ClusterTopology_
56
57//
b69620f8 58}
59
60//______________________________________________________________________________
61void AliITSUClusterizer::SetSegmentation(const AliITSUSegmentationPix *segm)
62{
63 // attach segmentation, if needed, reinitialize array
64 fSegm = segm;
65 SetNRow(fSegm->GetNRow()); // reinitialize if needed
66
67}
68
69//______________________________________________________________________________
70void AliITSUClusterizer::SetNRow(Int_t nr)
71{
72 // update buffers
73 int nrOld = GetUniqueID();
74 if (nrOld>=nr) return;
75 SetUniqueID(nr);
76 while (fDigitFreelistBptrFirst) {
77 AliITSUClusterizerClusterDigit *next = fDigitFreelistBptrFirst[kDigitChunkSize-1].fNext;
78 delete[] fDigitFreelistBptrFirst;
79 fDigitFreelistBptrFirst=next;
80 }
81 delete[] fPartFreelistBptr;
82 delete[] fCandFreelistBptr;
83 //
84 fPartFreelist=fPartFreelistBptr = new AliITSUClusterizerClusterPart[nr+1];
85 fCandFreelist=fCandFreelistBptr = new AliITSUClusterizerClusterCand[nr+1];
86 for (int i=0;i<nr;++i) {
87 fPartFreelistBptr[i].fNextInRow = &fPartFreelistBptr[i+1];
88 fCandFreelistBptr[i].fNext = &fCandFreelistBptr[i+1];
89 }
90}
91
92//______________________________________________________________________________
93AliITSUClusterizer::~AliITSUClusterizer()
94{
95 // d-tor
96 while (fDigitFreelistBptrFirst) {
97 AliITSUClusterizerClusterDigit *next = fDigitFreelistBptrFirst[kDigitChunkSize-1].fNext;
98 delete[] fDigitFreelistBptrFirst;
99 fDigitFreelistBptrFirst=next;
100 }
101 delete[] fPartFreelistBptr;
102 delete[] fCandFreelistBptr;
6cae87c5 103 //
b69620f8 104}
105
106//______________________________________________________________________________
107AliITSUClusterizer::AliITSUClusterizerClusterDigit* AliITSUClusterizer::AllocDigitFreelist()
108{
109 // allocate aux space
110 AliITSUClusterizerClusterDigit *tmp = new AliITSUClusterizerClusterDigit[kDigitChunkSize];
111 for (int i=0;i<kDigitChunkSize-2;++i) tmp[i].fNext=&tmp[i+1];
112 tmp[kDigitChunkSize-2].fNext=0;
113 tmp[kDigitChunkSize-1].fNext=0;
114 if (!fDigitFreelistBptrFirst) fDigitFreelistBptrFirst=tmp;
115 else fDigitFreelistBptrLast[kDigitChunkSize-1].fNext=tmp;
116 fDigitFreelistBptrLast=tmp;
117 return tmp;
118}
119
120//______________________________________________________________________________
121AliITSUClusterizer::AliITSUClusterizerClusterDigit* AliITSUClusterizer::NextDigit()
122{
123 // get next digit
c7811d12 124 while (fInputDigitsReadIndex<fInputDigits->GetEntriesFast()) {
125 AliITSUDigitPix *tmp=static_cast<AliITSUDigitPix*>(fInputDigits->UncheckedAt(fInputDigitsReadIndex++));
126 if (TMath::Abs(tmp->GetROCycle())>=fRecoParam->GetMaxROCycle()) continue;
b69620f8 127 AliITSUClusterizerClusterDigit *digit=AllocDigit();
128 digit->fDigit=tmp;
b69620f8 129 return digit;
130 }
c7811d12 131 return 0;
b69620f8 132}
133
134//______________________________________________________________________________
135void AliITSUClusterizer::AttachPartToCand(AliITSUClusterizerClusterCand *cand,AliITSUClusterizerClusterPart *part)
136{
137 // attach part
138 part->fParent = cand;
139 part->fPrevInCluster = 0;
140 part->fNextInCluster = cand->fFirstPart;
141 if (cand->fFirstPart) cand->fFirstPart->fPrevInCluster = part;
142 cand->fFirstPart=part;
143}
144
145//______________________________________________________________________________
146void AliITSUClusterizer::MergeCands(AliITSUClusterizerClusterCand *a,AliITSUClusterizerClusterCand *b)
147{
148 // merge cluster parts
149 AliITSUClusterizerClusterPart *ipart=b->fFirstPart;
150 AliITSUClusterizerClusterPart *jpart;
151 do {
152 jpart=ipart;
153 jpart->fParent=a;
154 } while ((ipart=ipart->fNextInCluster));
155 jpart->fNextInCluster=a->fFirstPart;
156 jpart->fNextInCluster->fPrevInCluster=jpart;
157 a->fFirstPart=b->fFirstPart;
158 // merge digits
159 b->fLastDigit->fNext=a->fFirstDigit;
160 a->fFirstDigit=b->fFirstDigit;
161 // DeallocCand(b);
162}
163
164//______________________________________________________________________________
5e375bb4 165void AliITSUClusterizer::Transform(AliITSUClusterPix *cluster,AliITSUClusterizerClusterCand *cand)
b69620f8 166{
5e375bb4 167 // convert set of digits to cluster data in LOCAL frame
168 const double k1to12 = 1./12;
4bac12be 169 static int maxLbinDigit = AliITSdigit::GetNTracks();
5e375bb4 170 //
4bac12be 171 fNLabels = 0;
b69620f8 172 Int_t n=0;
173 cand->fLastDigit->fNext=0;
5e375bb4 174 double x=0,z=0,xmn=1e9,xmx=-1e9,zmn=1e9,zmx=-1e9,px=0,pz=0;
175 float cx,cz;
29ad4146 176 int charge=0;
6cae87c5 177 //
b69620f8 178 for (AliITSUClusterizerClusterDigit *idigit=cand->fFirstDigit;idigit;idigit=idigit->fNext) {
4bac12be 179 AliITSdigit* dig = idigit->fDigit;
ee58ce21 180 fSegm->DetToLocal(dig->GetCoord2(),dig->GetCoord1(),cx,cz); // center of pixel
181 //
182 // account for possible diod shift
183 double ddx,ddz, dx=fSegm->Dpx(dig->GetCoord2()), dz=fSegm->Dpz(dig->GetCoord1());
184 fSegm->GetDiodShift(dig->GetCoord2(),dig->GetCoord1(),ddx,ddz);
185 //
29ad4146 186 charge += dig->GetSignal();
5e375bb4 187 x += cx;
188 z += cz;
189 if (cx<xmn) xmn=cx;
190 if (cx>xmx) xmx=cx;
191 if (cz<zmn) zmn=cz;
192 if (cz>zmx) zmx=cz;
ee58ce21 193 x += ddx*dx;
194 z += ddz*dz;
195 px += dx;
196 pz += dz;
4bac12be 197 //
198 if (!fRawData) {
199 for(Int_t dlab=0;dlab<maxLbinDigit;dlab++){
200 Int_t digitlab = (dig->GetTracks())[dlab];
201 if(digitlab<0) continue;
202 AddLabel(digitlab);
203 }
204 }
205 //
b69620f8 206 ++n;
207 }
b7e5fec5 208 Int_t nx=1,nz=1;
5e375bb4 209 double dx = xmx-xmn, dz = zmx-zmn;
210 if (n>1) {
211 double fac=1./n;
212 x *= fac; // mean coordinates
213 z *= fac;
214 px *= fac; // mean pitch
215 pz *= fac;
216 nx = 1+Nint(dx/px);
217 nz = 1+Nint(dz/pz);
b69620f8 218 }
889b1493 219 x -= fLorAngCorrection; // LorentzAngle correction
b69620f8 220 cluster->SetX(x);
b7e5fec5 221 cluster->SetZ(z);
5e375bb4 222 cluster->SetY(0);
c8d1f258 223 cluster->SetSigmaZ2(nz>1 ? dz*dz*k1to12 : pz*pz*k1to12);
224 cluster->SetSigmaY2(nx>1 ? dx*dx*k1to12 : px*px*k1to12);
5e375bb4 225 cluster->SetSigmaYZ(0);
226 cluster->SetFrameLoc();
b7e5fec5 227 if (n>AliITSUClusterPix::kMaskNPix) n = AliITSUClusterPix::kMaskNPix;
228 cluster->SetNxNzN(nx>AliITSUClusterPix::kMaskNX ? AliITSUClusterPix::kMaskNX : nx,
229 nz>AliITSUClusterPix::kMaskNZ ? AliITSUClusterPix::kMaskNZ : nz,
230 n>AliITSUClusterPix::kMaskNPix ? AliITSUClusterPix::kMaskNPix : n);
29ad4146 231 cluster->SetQ(charge); // note: this is MC info
b7e5fec5 232 // if (cluster->GetNPix()!=n || cluster->GetNx()!=nx || cluster->GetNz()!=nz) {
233 // printf("PROBLEM: %d %d %d -> %d %d %d\n",nx,nz,n,
234 // cluster->GetNx(),cluster->GetNz(),cluster->GetNPix());
235 // cluster->Dump();
236 // }
5e375bb4 237 //
4bac12be 238 if (!fRawData) {
239 CheckLabels();
240 int nl = Min(kMaxLabInCluster,fNLabels);
241 for (int i=nl;i--;) cluster->SetLabel(fCurrLabels[i],i);
242 }
243 //
b69620f8 244 // Set Volume id
245 cluster->SetVolumeId(fVolID);
246 // printf("mod %d: (%.4lf,%.4lf)cm\n",fVolID,x,z);
6cae87c5 247 //
248#ifdef _ClusterTopology_
6f54924d 249 FillClusterTopology(cand,cluster);
6cae87c5 250#endif //_ClusterTopology_
251 //
b69620f8 252}
253
254//______________________________________________________________________________
255void AliITSUClusterizer::CloseCand(AliITSUClusterizerClusterCand *cand)
256{
257 // finish cluster
5e375bb4 258 AliITSUClusterPix *cluster = (AliITSUClusterPix*)NextCluster();
b69620f8 259 Transform(cluster,cand);
260 DeallocDigits(cand->fFirstDigit,cand->fLastDigit);
261 DeallocCand(cand);
262}
263
264//______________________________________________________________________________
265void AliITSUClusterizer::ClosePart(AliITSUClusterizerClusterPart *part)
266{
267 // finish cluster part
268 AliITSUClusterizerClusterCand *cand=part->fParent;
269 DetachPartFromCand(cand,part);
270 DeallocPart(part);
271 if (!cand->fFirstPart) CloseCand(cand);
272}
273
274//______________________________________________________________________________
275void AliITSUClusterizer::CloseRemainingParts(AliITSUClusterizerClusterPart *part)
276{
277 // finish what is left
278 while (part) {
279 AliITSUClusterizerClusterPart *next=part->fNextInRow;
280 ClosePart(part);
281 part=next;
282 }
283}
284
285//______________________________________________________________________________
286void AliITSUClusterizer::Clusterize()
287{
4bac12be 288 // main algo for MC clustererization
289 SetRawData(kFALSE);
290 //
b69620f8 291 AliITSUClusterizerClusterDigit *iDigit=NextDigit();
292 AliITSUClusterizerClusterPart *iPrevRowBegin=0;
293 AliITSUClusterizerClusterPart *iNextRowBegin=0;
294 AliITSUClusterizerClusterPart *iPrevRow=0;
295 AliITSUClusterizerClusterPart *iNextRow=0;
296 Int_t lastV=0;
6cae87c5 297 //
b69620f8 298 while (iDigit) {
4bac12be 299 if (iDigit->fDigit->GetCoord2()!=lastV) {
b69620f8 300 // NEW ROW
301 if (iNextRow) iNextRow->fNextInRow=0;
302 if (iPrevRowBegin) CloseRemainingParts(iPrevRowBegin);
4bac12be 303 if (iDigit->fDigit->GetCoord2()==lastV+1) {
b69620f8 304 iPrevRowBegin=iNextRowBegin;
305 iPrevRow =iNextRowBegin;
306 }
307 else {
308 // there was an empty row
309 CloseRemainingParts(iNextRowBegin);
310 iPrevRowBegin=0;
311 iPrevRow =0;
312 }
313 iNextRowBegin=0;
314 iNextRow =0;
4bac12be 315 lastV=iDigit->fDigit->GetCoord2();
b69620f8 316 }
317 // skip cluster parts before this digit
ee58ce21 318 int limCol = iDigit->fDigit->GetCoord1()-fAllowDiagonalClusterization;
319 while (iPrevRow && iPrevRow->fUEnd < limCol) {
b69620f8 320 iPrevRow=iPrevRow->fNextInRow;
321 }
322 // find the longest continous line of digits [iDigit,pDigit]=[iDigit,jDigit)
323 AliITSUClusterizerClusterCand *cand=AllocCand();
324 AliITSUClusterizerClusterDigit *pDigit=iDigit;
325 AliITSUClusterizerClusterDigit *jDigit=NextDigit();
326 cand->fFirstPart=0;
327 cand->fFirstDigit=cand->fLastDigit=iDigit; // NB: first diggit is attached differently
328 iDigit->fNext=0;
4bac12be 329 Int_t lastU =iDigit->fDigit->GetCoord1();
b69620f8 330 Int_t lastU1=lastU+1;
4bac12be 331 while (jDigit && jDigit->fDigit->GetCoord1()==lastU1 && jDigit->fDigit->GetCoord2()==lastV) {
b69620f8 332 pDigit=jDigit;
333 jDigit=NextDigit();
334 AttachDigitToCand(cand,pDigit);
335 ++lastU1;
336 }
ee58ce21 337 if (!fAllowDiagonalClusterization) --lastU1;
b69620f8 338 AliITSUClusterizerClusterPart *part=AllocPart();
339 part->fUBegin=lastU ;
340 part->fUEnd =lastU1;
341 AttachPartToCand(cand,part);
342 // merge all cluster candidates of the previous line touching this one,
343 // advance to the last one, but keep that one the next active one
344 AliITSUClusterizerClusterPart *jPrevRow=iPrevRow;
345 while (jPrevRow && jPrevRow->fUBegin<=lastU1) {
346 if (jPrevRow->fParent!=cand) {
347 MergeCands(jPrevRow->fParent,cand);
348 DeallocCand(cand);
349 cand=jPrevRow->fParent;
350 }
351 iPrevRow=jPrevRow;
352 jPrevRow=jPrevRow->fNextInRow;
353 }
354 if (iNextRow)
355 iNextRow->fNextInRow=part;
356 else
357 iNextRowBegin=part;
358 iNextRow=part;
359 iDigit=jDigit;
360 }
361 // remove remaining cluster parts
362 CloseRemainingParts(iPrevRowBegin);
363 if (iNextRow) iNextRow->fNextInRow=0;
364 CloseRemainingParts(iNextRowBegin);
365 return;
366}
889b1493 367
368//______________________________________________________________________________
369void AliITSUClusterizer::PrepareLorentzAngleCorrection(Double_t bz)
370{
371 // calculate parameters for Lorentz Angle correction. Must be called
372 // after setting segmentation and recoparams
373 fLorAngCorrection = 0.5*fRecoParam->GetTanLorentzAngle(fLayerID)*bz/kNominalBz*fSegm->Dy();
374}
4bac12be 375
376//______________________________________________________________________
377void AliITSUClusterizer::CheckLabels()
378{
379 // Tries to find mother's labels
380 //
381 if (fNLabels<1) return;
382 AliRunLoader *rl = AliRunLoader::Instance();
383 if(!rl) return;
384 TTree *trK=(TTree*)rl->TreeK();
385 if (!trK) return;
386 //
387 static int labS[kMaxLabels];
388 static float kine[kMaxLabels];
389 Int_t nlabels = fNLabels;
390 Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
391 for (Int_t i=fNLabels;i--;) labS[i] = fCurrLabels[i];
392 //
393 for (Int_t i=0;i<nlabels;i++) {
394 Int_t label = labS[i];
395 if (label>=ntracks) continue;
396 TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
397 kine[i] = part->Energy() - part->GetCalcMass(); // kinetic energy
81a0c667 398 if (kine[i] < 0.001) { // reduce soft particles from the same cluster
4bac12be 399 Int_t m=part->GetFirstMother();
400 if (m<0) continue; // primary
401 //
402 if (part->GetStatusCode()>0) continue;
403 //
404 // if the parent is within the same cluster, assign parent's label
405 for (int j=0;j<nlabels;j++) if (labS[j]==m) { labS[i] = m; break;}
406 }
407 }
408 //
409 if (nlabels>kMaxLabInCluster) { // only 3 labels are stored in cluster, sort in decreasing momentum
410 static int ind[kMaxLabels],labSS[kMaxLabels];
411 TMath::Sort(nlabels,kine,ind);
412 for (int i=nlabels;i--;) labSS[i] = labS[i];
413 for (int i=nlabels;i--;) labS[i] = labSS[ind[i]];
414 }
415 //
416 //compress labels -- if multi-times the same
417 for (Int_t i=0;i<nlabels;i++) fCurrLabels[i]=-2;
418 fNLabels = 0;
419 int j=0;
420 for (int i=0;i<nlabels;i++) {
421 for (j=fNLabels;j--;) if (labS[i]==fCurrLabels[j]) break; // the label already there
422 if (j<0) fCurrLabels[fNLabels++] = labS[i];
423 }
424 //
425}
6cae87c5 426
427//------------------------------------------------------------------------
428
429#ifdef _ClusterTopology_
430//
431//______________________________________________________________________
6f54924d 432void AliITSUClusterizer::FillClusterTopology(const AliITSUClusterizerClusterCand *cand, AliITSUClusterPix* cl) const
6cae87c5 433{
6f54924d 434 // fill cluster topology bit pattern
6cae87c5 435 //
6cae87c5 436 //
6f54924d 437 UShort_t mnCol=0xffff,mxCol=0,mnRow=0xffff,mxRow=0;
6cae87c5 438 for (AliITSUClusterizerClusterDigit *idigit=cand->fFirstDigit;idigit;idigit=idigit->fNext) {
439 AliITSdigit* dig = idigit->fDigit;
6f54924d 440 if (mnCol>dig->GetCoord1()) mnCol = (UShort_t)dig->GetCoord1();
441 if (mxCol<dig->GetCoord1()) mxCol = (UShort_t)dig->GetCoord1();
442 if (mnRow>dig->GetCoord2()) mnRow = (UShort_t)dig->GetCoord2();
443 if (mxRow<dig->GetCoord2()) mxRow = (UShort_t)dig->GetCoord2();
6cae87c5 444 }
6f54924d 445 UShort_t nColSpan = UShort_t(mxCol-mnCol+1);
446 UShort_t nRowSpan = UShort_t(mxRow-mnRow+1);
447 UShort_t nColSpanW = nColSpan;
448 UShort_t nRowSpanW = nRowSpan;
449 if (nColSpan*nRowSpan>AliITSUClusterPix::kMaxPatternBits) { // need to store partial info
450 // will crtail largest dimension
451 if (nColSpan>nRowSpan) {
452 if ( (nColSpanW=AliITSUClusterPix::kMaxPatternBits/nRowSpan)==0 ) {
453 nColSpanW = 1;
454 nRowSpanW = AliITSUClusterPix::kMaxPatternBits;
455 }
456 }
457 else {
458 if ( (nRowSpanW=AliITSUClusterPix::kMaxPatternBits/nColSpan)==0 ) {
459 nRowSpanW = 1;
460 nColSpanW = AliITSUClusterPix::kMaxPatternBits;
461 }
6cae87c5 462 }
6f54924d 463 }
464 cl->ResetPattern();
465 //
466 cl->SetPatternRowSpan(nRowSpanW,nRowSpanW<nRowSpan);
467 cl->SetPatternColSpan(nColSpanW,nColSpanW<nColSpan);
468 cl->SetPatternMinRow(mnRow);
469 cl->SetPatternMinCol(mnCol);
470 //
471 for (AliITSUClusterizerClusterDigit *idigit=cand->fFirstDigit;idigit;idigit=idigit->fNext) {
472 AliITSdigit* dig = idigit->fDigit;
473 UInt_t ir = dig->GetCoord2()-mnRow;
474 UInt_t ic = dig->GetCoord1()-mnCol;
475 if (ir<nRowSpanW && ic<nColSpanW) cl->SetPixel(ir,ic);
6cae87c5 476 }
477 //
6cae87c5 478}
479
480#endif //_ClusterTopology_