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