Classes moved to STEERBase.
[u/mrichter/AliRoot.git] / TRD / AliTRDseed.cxx
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  **************************************************************************/
15                                                       
16 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 //                                                                           //
20 //  The seed of a local TRD track                                            //  
21 //                                                                           //
22 ///////////////////////////////////////////////////////////////////////////////
23
24 #include "TMath.h"
25 #include "TLinearFitter.h"
26
27 #include "AliMathBase.h"
28
29 #include "AliTRDseed.h"
30 #include "AliTRDcluster.h"
31 #include "AliTRDtracker.h"
32
33 ClassImp(AliTRDseed)
34
35 //_____________________________________________________________________________
36 AliTRDseed::AliTRDseed() 
37   :TObject()
38   ,fTilt(0)
39   ,fPadLength(0)
40   ,fX0(0)
41   ,fSigmaY(0)
42   ,fSigmaY2(0)
43   ,fMeanz(0)
44   ,fZProb(0)
45   ,fN(0)
46   ,fN2(0)
47   ,fNUsed(0)
48   ,fFreq(0)
49   ,fNChange(0)
50   ,fMPads(0)
51   ,fC(0)
52   ,fCC(0)
53   ,fChi2(1.0e10)
54   ,fChi2Z(1.0e10)
55 {
56   //
57   // Default constructor  
58   //
59
60   for (Int_t i = 0; i < 25; i++) {
61     fX[i]        = 0;   // x position
62     fY[i]        = 0;   // y position
63     fZ[i]        = 0;   // z position
64     fIndexes[i]  = 0;   // Indexes
65     fClusters[i] = 0x0; // Clusters
66     fUsable[i]   = 0;   // Indication  - usable cluster
67   }
68
69   for (Int_t i = 0; i < 2; i++) {
70     fYref[i]     = 0;   // Reference y
71     fZref[i]     = 0;   // Reference z
72     fYfit[i]     = 0;   // Y fit position +derivation
73     fYfitR[i]    = 0;   // Y fit position +derivation
74     fZfit[i]     = 0;   // Z fit position
75     fZfitR[i]    = 0;   // Z fit position
76     fLabels[i]   = 0;   // Labels
77   }
78
79 }
80
81 //_____________________________________________________________________________
82 AliTRDseed::AliTRDseed(const AliTRDseed &s)
83   :TObject(s)
84   ,fTilt(s.fTilt)
85   ,fPadLength(s.fPadLength)
86   ,fX0(s.fX0)
87   ,fSigmaY(s.fSigmaY)
88   ,fSigmaY2(s.fSigmaY2)
89   ,fMeanz(s.fMeanz)
90   ,fZProb(s.fZProb)
91   ,fN(s.fN)
92   ,fN2(s.fN2)
93   ,fNUsed(s.fNUsed)
94   ,fFreq(s.fFreq)
95   ,fNChange(s.fNChange)
96   ,fMPads(s.fMPads)
97   ,fC(s.fC)
98   ,fCC(s.fCC)
99   ,fChi2(s.fChi2)
100   ,fChi2Z(s.fChi2Z)
101 {
102   //
103   // Copy constructor  
104   //
105
106   for (Int_t i = 0; i < 25; i++) {
107     fX[i]        = s.fX[i];        // x position
108     fY[i]        = s.fY[i];        // y position
109     fZ[i]        = s.fZ[i];        // z position
110     fIndexes[i]  = s.fIndexes[i];  // Indexes
111     fClusters[i] = s.fClusters[i]; // Clusters
112     fUsable[i]   = s.fUsable[i];   // Indication  - usable cluster
113   }
114
115   for (Int_t i = 0; i < 2; i++) {
116     fYref[i]     = s.fYref[i];     // Reference y
117     fZref[i]     = s.fZref[i];     // Reference z
118     fYfit[i]     = s.fYfit[i];     // Y fit position +derivation
119     fYfitR[i]    = s.fYfitR[i];    // Y fit position +derivation
120     fZfit[i]     = s.fZfit[i];     // Z fit position
121     fZfitR[i]    = s.fZfitR[i];    // Z fit position
122     fLabels[i]   = s.fLabels[i];   // Labels
123   }
124
125 }
126
127 //_____________________________________________________________________________
128 void AliTRDseed::Reset()
129 {
130   //
131   // Reset seed
132   //
133
134   for (Int_t i = 0; i < 25; i++) {
135     fX[i]        = 0;  // X position
136     fY[i]        = 0;  // Y position
137     fZ[i]        = 0;  // Z position
138     fIndexes[i]  = 0;  // Indexes
139     fClusters[i] = 0;  // Clusters
140     fUsable[i]   = kFALSE;    
141   }
142
143   for (Int_t i = 0; i < 2; i++) {
144     fYref[i]     = 0;  // Reference y
145     fZref[i]     = 0;  // Reference z
146     fYfit[i]     = 0;  // Y fit position +derivation
147     fYfitR[i]    = 0;  // Y fit position +derivation
148     fZfit[i]     = 0;  // Z fit position
149     fZfitR[i]    = 0;  // Z fit position
150     fLabels[i]   = -1; // Labels
151   }
152   fSigmaY  = 0;        // "Robust" sigma in y
153   fSigmaY2 = 0;        // "Robust" sigma in y
154   fMeanz   = 0;        // Mean vaue of z
155   fZProb   = 0;        // Max probbable z
156   fMPads   = 0;
157   fN       = 0;        // Number of associated clusters
158   fN2      = 0;        // Number of not crossed
159   fNUsed   = 0;        // Number of used clusters
160   fNChange = 0;        // Change z counter
161
162 }
163
164 //_____________________________________________________________________________
165 void AliTRDseed::CookLabels()
166 {
167   //
168   // Cook 2 labels for seed
169   //
170
171   Int_t labels[200];
172   Int_t out[200];
173   Int_t nlab = 0;
174
175   for (Int_t i = 0; i < 25; i++) {
176     if (!fClusters[i]) continue;
177     for (Int_t ilab = 0; ilab < 3; ilab++) {
178       if (fClusters[i]->GetLabel(ilab) >= 0) {
179         labels[nlab] = fClusters[i]->GetLabel(ilab);
180         nlab++;
181       }
182     }
183   }
184
185   Int_t nlab2 = AliTRDtracker::Freq(nlab,labels,out,kTRUE);
186   fLabels[0] = out[0];
187   if ((nlab2  > 1) && 
188       (out[3] > 1)) {
189     fLabels[1] = out[2];
190   }
191
192 }
193
194 //_____________________________________________________________________________
195 void AliTRDseed::UseClusters()
196 {
197   //
198   // Use clusters
199   //
200
201   for (Int_t i = 0; i < 25; i++) {
202     if (!fClusters[i]) continue;
203     if (!(fClusters[i]->IsUsed())) fClusters[i]->Use();
204   }
205
206 }
207
208 //_____________________________________________________________________________
209 void AliTRDseed::Update()
210 {
211   //
212   // Update the seed.
213   //
214
215   const Float_t kRatio  = 0.8;
216   const Int_t   kClmin  = 6;
217   const Float_t kmaxtan = 2;
218
219   if (TMath::Abs(fYref[1]) > kmaxtan) return;              // Track inclined too much
220
221   Float_t  sigmaexp  = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction
222   Float_t  ycrosscor = fPadLength * fTilt * 0.5;           // Y correction for crossing 
223   fNChange = 0;
224
225   Double_t sumw;
226   Double_t sumwx;
227   Double_t sumwx2;
228   Double_t sumwy;
229   Double_t sumwxy;
230   Double_t sumwz;
231   Double_t sumwxz;
232
233   Int_t    zints[25];                    // Histograming of the z coordinate 
234                                          // Get 1 and second max probable coodinates in z
235   Int_t    zouts[50];       
236   Float_t  allowedz[25];                 // Allowed z for given time bin
237   Float_t  yres[25];                     // Residuals from reference
238   Float_t  anglecor = fTilt * fZref[1];  // Correction to the angle
239   
240   
241   fN  = 0; 
242   fN2 = 0;
243   for (Int_t i = 0; i < 25; i++) {
244     yres[i] = 10000.0;
245     if (!fClusters[i]) continue;
246     yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i];   // Residual y
247     zints[fN] = Int_t(fZ[i]);
248     fN++;    
249   }
250
251   if (fN < kClmin) return;
252   Int_t nz = AliTRDtracker::Freq(fN,zints,zouts,kFALSE);
253   fZProb   = zouts[0];
254   if (nz <= 1) zouts[3] = 0;
255   if (zouts[1] + zouts[3] < kClmin) return;
256   
257   // Z distance bigger than pad - length
258   if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) {
259     zouts[3]=0;           
260   }
261   
262   Int_t  breaktime = -1;
263   Bool_t mbefore   = kFALSE;
264   Int_t  cumul[25][2];
265   Int_t  counts[2] = { 0, 0 };
266   
267   if (zouts[3] >= 3) {
268
269     //
270     // Find the break time allowing one chage on pad-rows
271     // with maximal numebr of accepted clusters
272     //
273     fNChange = 1;
274     for (Int_t i = 0; i < 25; i++) {
275       cumul[i][0] = counts[0];
276       cumul[i][1] = counts[1];
277       if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++;
278       if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++;
279     }
280     Int_t  maxcount = 0;
281     for (Int_t i = 0; i < 24; i++) {
282       Int_t after  = cumul[24][0] - cumul[i][0];
283       Int_t before = cumul[i][1];
284       if (after + before > maxcount) { 
285         maxcount  = after + before; 
286         breaktime = i;
287         mbefore   = kFALSE;
288       }
289       after  = cumul[24][1] - cumul[i][1];
290       before = cumul[i][0];
291       if (after + before > maxcount) { 
292         maxcount  = after + before; 
293         breaktime = i;
294         mbefore   = kTRUE;
295       }
296     }
297
298     breaktime -= 1;
299
300   }
301
302   for (Int_t i = 0; i < 25; i++) {
303     if (i >  breaktime) allowedz[i] =   mbefore  ? zouts[2] : zouts[0];
304     if (i <= breaktime) allowedz[i] = (!mbefore) ? zouts[2] : zouts[0];
305   }  
306
307   if (((allowedz[0] > allowedz[24]) && (fZref[1] < 0)) || 
308       ((allowedz[0] < allowedz[24]) && (fZref[1] > 0))) {
309     //
310     // Tracklet z-direction not in correspondance with track z direction 
311     //
312     fNChange = 0;
313     for (Int_t i = 0; i < 25; i++) {
314       allowedz[i] = zouts[0];  // Only longest taken
315     } 
316   }
317   
318   if (fNChange > 0) {
319     //
320     // Cross pad -row tracklet  - take the step change into account
321     //
322     for (Int_t i = 0; i < 25; i++) {
323       if (!fClusters[i]) continue; 
324       if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
325       yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i];   // Residual y
326       if (TMath::Abs(fZ[i] - fZProb) > 2) {
327         if (fZ[i] > fZProb) yres[i] += fTilt * fPadLength;
328         if (fZ[i] < fZProb) yres[i] -= fTilt * fPadLength;
329       }
330     }
331   }
332   
333   Double_t yres2[25];
334   Double_t mean;
335   Double_t sigma;
336   for (Int_t i = 0; i < 25; i++) {
337     if (!fClusters[i]) continue;
338     if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
339     yres2[fN2] = yres[i];
340     fN2++;
341   }
342   if (fN2 < kClmin) {
343     fN2 = 0;
344     return;
345   }
346   AliMathBase::EvaluateUni(fN2,yres2,mean,sigma,Int_t(fN2*kRatio-2));
347   if (sigma < sigmaexp * 0.8) {
348     sigma = sigmaexp;
349   }
350   fSigmaY = sigma;
351
352   // Reset sums
353   sumw   = 0; 
354   sumwx  = 0; 
355   sumwx2 = 0;
356   sumwy  = 0; 
357   sumwxy = 0; 
358   sumwz  = 0;
359   sumwxz = 0;
360
361   fN2    = 0;
362   fMeanz = 0;
363   fMPads = 0;
364
365   for (Int_t i = 0; i < 25; i++) {
366
367     fUsable[i] = kFALSE;
368     if (!fClusters[i]) continue;
369     if (TMath::Abs(fZ[i] - allowedz[i]) > 2)  continue;
370     if (TMath::Abs(yres[i] - mean) > 4.0 * sigma) continue;
371     fUsable[i] = kTRUE;
372     fN2++;
373     fMPads += fClusters[i]->GetNPads();
374     Float_t weight = 1.0;
375     if (fClusters[i]->GetNPads() > 4) weight = 0.5;
376     if (fClusters[i]->GetNPads() > 5) weight = 0.2;
377    
378     Double_t x = fX[i];
379     sumw   += weight; 
380     sumwx  += x * weight; 
381     sumwx2 += x*x * weight;
382     sumwy  += weight * yres[i];  
383     sumwxy += weight * (yres[i]) * x;
384     sumwz  += weight * fZ[i];    
385     sumwxz += weight * fZ[i] * x;
386
387   }
388
389   if (fN2 < kClmin){
390     fN2 = 0;
391     return;
392   }
393   fMeanz = sumwz / sumw;
394   Float_t correction = 0;
395   if (fNChange > 0) {
396     // Tracklet on boundary
397     if (fMeanz < fZProb) correction =  ycrosscor;
398     if (fMeanz > fZProb) correction = -ycrosscor;
399   }
400
401   Double_t det = sumw * sumwx2 - sumwx * sumwx;
402   fYfitR[0]    = (sumwx2 * sumwy  - sumwx * sumwxy) / det;
403   fYfitR[1]    = (sumw   * sumwxy - sumwx * sumwy)  / det;
404   
405   fSigmaY2 = 0;
406   for (Int_t i = 0; i < 25; i++) {    
407     if (!fUsable[i]) continue;
408     Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i];
409     fSigmaY2 += delta*delta;
410   }
411   fSigmaY2 = TMath::Sqrt(fSigmaY2 / Float_t(fN2-2));
412   
413   fZfitR[0]  = (sumwx2 * sumwz  - sumwx * sumwxz) / det;
414   fZfitR[1]  = (sumw   * sumwxz - sumwx * sumwz)  / det;
415   fZfit[0]   = (sumwx2 * sumwz  - sumwx * sumwxz) / det;
416   fZfit[1]   = (sumw   * sumwxz - sumwx * sumwz)  / det;
417   fYfitR[0] += fYref[0] + correction;
418   fYfitR[1] += fYref[1];
419   fYfit[0]   = fYfitR[0];
420   fYfit[1]   = fYfitR[1];
421     
422   UpdateUsed();
423
424 }
425
426 //_____________________________________________________________________________
427 void AliTRDseed::UpdateUsed()
428 {
429   //
430   // Update used seed
431   //
432
433   fNUsed = 0;
434   for (Int_t i = 0; i < 25; i++) {
435     if (!fClusters[i]) {
436       continue;
437     }
438     if ((fClusters[i]->IsUsed())) {
439       fNUsed++;
440     }
441   }
442
443 }
444
445 //_____________________________________________________________________________
446 Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror)
447 {
448   //
449   // Fit the Rieman tilt
450   //
451
452   // Fitting with tilting pads - kz not fixed
453   TLinearFitter fitterT2(4,"hyp4");  
454   fitterT2.StoreData(kTRUE);
455   Float_t xref2 = (cseed[2].fX0 + cseed[3].fX0) * 0.5; // Reference x0 for z
456   
457   Int_t npointsT = 0;
458   fitterT2.ClearPoints();
459
460   for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
461
462     if (!cseed[iLayer].IsOK()) continue;
463     Double_t tilt = cseed[iLayer].fTilt;
464
465     for (Int_t itime = 0; itime < 25; itime++) {
466
467       if (!cseed[iLayer].fUsable[itime]) continue;
468       // x relative to the midle chamber
469       Double_t x = cseed[iLayer].fX[itime] + cseed[iLayer].fX0 - xref2;  
470       Double_t y = cseed[iLayer].fY[itime];
471       Double_t z = cseed[iLayer].fZ[itime];
472
473       //
474       // Tilted rieman
475       //
476       Double_t uvt[6];
477       Double_t x2 = cseed[iLayer].fX[itime] + cseed[iLayer].fX0;      // Global x
478       Double_t t  = 1.0 / (x2*x2 + y*y);
479       uvt[1]  = t;
480       uvt[0]  = 2.0 * x2   * uvt[1];
481       uvt[2]  = 2.0 * tilt * uvt[1];
482       uvt[3]  = 2.0 * tilt *uvt[1] * x;       
483       uvt[4]  = 2.0 * (y + tilt * z) * uvt[1];
484       
485       Double_t error = 2.0 * uvt[1];
486       if (terror) {
487         error *= cseed[iLayer].fSigmaY;
488       }
489       else {
490         error *= 0.2; //Default error
491       } 
492       fitterT2.AddPoint(uvt,uvt[4],error);
493       npointsT++;
494
495     }
496
497   }
498
499   fitterT2.Eval();
500   Double_t rpolz0 = fitterT2.GetParameter(3);
501   Double_t rpolz1 = fitterT2.GetParameter(4);       
502
503   //
504   // Linear fitter  - not possible to make boundaries
505   // non accept non possible z and dzdx combination
506   //        
507   Bool_t acceptablez = kTRUE;
508   for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
509     if (cseed[iLayer].IsOK()) {
510       Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].fX0 - xref2);
511       if (TMath::Abs(cseed[iLayer].fZProb - zT2) > cseed[iLayer].fPadLength * 0.5 + 1.0) {
512         acceptablez = kFALSE;
513       }
514     }
515   }
516   if (!acceptablez) {
517     Double_t zmf  = cseed[2].fZref[0] + cseed[2].fZref[1] * (xref2 - cseed[2].fX0);
518     Double_t dzmf = (cseed[2].fZref[1] + cseed[3].fZref[1]) * 0.5;
519     fitterT2.FixParameter(3,zmf);
520     fitterT2.FixParameter(4,dzmf);
521     fitterT2.Eval();
522     fitterT2.ReleaseParameter(3);
523     fitterT2.ReleaseParameter(4);
524     rpolz0 = fitterT2.GetParameter(3);
525     rpolz1 = fitterT2.GetParameter(4);
526   }
527   
528   Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT);  
529   Double_t params[3];
530   params[0] =  fitterT2.GetParameter(0);
531   params[1] =  fitterT2.GetParameter(1);
532   params[2] =  fitterT2.GetParameter(2);            
533   Double_t curvature =  1.0 + params[1] * params[1] - params[2] * params[0];
534
535   for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
536
537     Double_t  x  = cseed[iLayer].fX0;
538     Double_t  y  = 0;
539     Double_t  dy = 0;
540     Double_t  z  = 0;
541     Double_t  dz = 0;
542
543     // y
544     Double_t res2 = (x * params[0] + params[1]);
545     res2 *= res2;
546     res2  = 1.0 - params[2]*params[0] + params[1]*params[1] - res2;
547     if (res2 >= 0) {
548       res2 = TMath::Sqrt(res2);
549       y    = (1.0 - res2) / params[0];
550     }
551
552     //dy
553     Double_t x0 = -params[1] / params[0];
554     if (-params[2]*params[0] + params[1]*params[1] + 1 > 0) {
555       Double_t rm1 = params[0] / TMath::Sqrt(-params[2]*params[0] + params[1]*params[1] + 1); 
556       if (1.0/(rm1*rm1) - (x-x0) * (x-x0) > 0.0) {
557         Double_t res = (x - x0) / TMath::Sqrt(1.0 / (rm1*rm1) - (x-x0)*(x-x0));
558         if (params[0] < 0) res *= -1.0;
559         dy = res;
560       }
561     }
562     z  = rpolz0 + rpolz1 * (x - xref2);
563     dz = rpolz1;
564     cseed[iLayer].fYref[0] = y;
565     cseed[iLayer].fYref[1] = dy;
566     cseed[iLayer].fZref[0] = z;
567     cseed[iLayer].fZref[1] = dz;
568     cseed[iLayer].fC       = curvature;
569     
570   }
571
572   return chi2TR;
573
574 }