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