Updates for the tracking code for calibration/alignment awareness
[u/mrichter/AliRoot.git] / TRD / AliTRDseedV1.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 TRD track seed                                                    //
21 //                                                                        //
22 //  Authors:                                                              //
23 //    Alex Bercuci <A.Bercuci@gsi.de>                                     //
24 //    Markus Fasel <M.Fasel@gsi.de>                                       //
25 //                                                                        //
26 ////////////////////////////////////////////////////////////////////////////
27
28 #include "TMath.h"
29 #include "TLinearFitter.h"
30
31 #include "AliLog.h"
32 #include "AliMathBase.h"
33
34 #include "AliTRDseedV1.h"
35 #include "AliTRDcluster.h"
36 #include "AliTRDtrack.h"
37 #include "AliTRDcalibDB.h"
38 #include "AliTRDstackLayer.h"
39 #include "AliTRDrecoParam.h"
40 #include "AliTRDgeometry.h"
41 #include "Cal/AliTRDCalPID.h"
42
43 #define SEED_DEBUG
44
45 ClassImp(AliTRDseedV1)
46
47 //____________________________________________________________________
48 AliTRDseedV1::AliTRDseedV1(Int_t layer, AliTRDrecoParam *p) 
49   :AliTRDseed()
50   ,fPlane(layer)
51   ,fOwner(kFALSE)
52   ,fMom(0.)
53   ,fSnp(0.)
54   ,fTgl(0.)
55   ,fdX(0.)
56   ,fRecoParam(p)
57 {
58   //
59   // Constructor
60   //
61         for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = 0.;
62         for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec]  = -1.;
63 }
64
65 //____________________________________________________________________
66 AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
67   :AliTRDseed((AliTRDseed&)ref)
68   ,fPlane(ref.fPlane)
69   ,fOwner(kFALSE)
70   ,fMom(ref.fMom)
71   ,fSnp(ref.fSnp)
72   ,fTgl(ref.fTgl)
73   ,fdX(ref.fdX)
74   ,fRecoParam(ref.fRecoParam)
75 {
76   //
77   // Copy Constructor performing a deep copy
78   //
79
80         //AliInfo("");
81         if(ref.fOwner) SetOwner();
82         for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = ref.fdEdx[islice];
83         for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = ref.fProb[ispec];
84 }
85
86
87 //____________________________________________________________________
88 AliTRDseedV1& AliTRDseedV1::operator=(const AliTRDseedV1 &ref)
89 {
90   //
91   // Assignment Operator using the copy function
92   //
93
94         //AliInfo("");
95         if(this != &ref){
96                 ref.Copy(*this);
97         }
98         return *this;
99
100 }
101
102 //____________________________________________________________________
103 AliTRDseedV1::~AliTRDseedV1()
104 {
105   //
106   // Destructor. The RecoParam object belongs to the underlying tracker.
107   //
108
109         //AliInfo(Form("fOwner[%s]", fOwner?"YES":"NO"));
110
111         if(fOwner) 
112                 for(int itb=0; itb<knTimebins; itb++){
113                         if(!fClusters[itb]) continue; 
114                         //AliInfo(Form("deleting c %p @ %d", fClusters[itb], itb));
115                         delete fClusters[itb];
116                         fClusters[itb] = 0x0;
117                 }
118 }
119
120 //____________________________________________________________________
121 void AliTRDseedV1::Copy(TObject &ref) const
122 {
123   //
124   // Copy function
125   //
126
127         //AliInfo("");
128         AliTRDseedV1 &target = (AliTRDseedV1 &)ref; 
129         
130         target.fPlane         = fPlane;
131         target.fMom           = fMom;
132         target.fSnp           = fSnp;
133         target.fTgl           = fTgl;
134         target.fdX            = fdX;
135         target.fRecoParam     = fRecoParam;
136         
137         for(int islice=0; islice < knSlices; islice++) target.fdEdx[islice] = fdEdx[islice];
138         for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) target.fProb[ispec] = fProb[ispec];
139         
140         AliTRDseed::Copy(target);
141 }
142
143
144 //____________________________________________________________
145 void AliTRDseedV1::Init(AliTRDtrack *track)
146 {
147 // Initialize this tracklet using the track information
148 //
149 // Parameters:
150 //   track - the TRD track used to initialize the tracklet
151 // 
152 // Detailed description
153 // The function sets the starting point and direction of the
154 // tracklet according to the information from the TRD track.
155 // 
156 // Caution
157 // The TRD track has to be propagated to the beginning of the
158 // chamber where the tracklet will be constructed
159 //
160
161         Double_t y, z; 
162         track->GetProlongation(fX0, y, z);
163         fYref[0] = y;
164         fYref[1] = track->GetSnp()/(1. - track->GetSnp()*track->GetSnp());
165         fZref[0] = z;
166         fZref[1] = track->GetTgl();
167
168         //printf("Tracklet ref x[%7.3f] y[%7.3f] z[%7.3f], snp[%f] tgl[%f]\n", fX0, fYref[0], fZref[0], track->GetSnp(), track->GetTgl());
169 }
170
171
172 //____________________________________________________________________
173 void AliTRDseedV1::CookdEdx(Int_t nslices)
174 {
175 // Calculates average dE/dx for all slices and store them in the internal array fdEdx. 
176 //
177 // Parameters:
178 //  nslices : number of slices for which dE/dx should be calculated
179 // Output:
180 //  store results in the internal array fdEdx. This can be accessed with the method
181 //  AliTRDseedV1::GetdEdx()
182 //
183 // Detailed description
184 // Calculates average dE/dx for all slices. Depending on the PID methode 
185 // the number of slices can be 3 (LQ) or 8(NN). 
186 // The calculation of dQ/dl are done using the tracklet fit results (see AliTRDseedV1::GetdQdl(Int_t)) i.e.
187 //
188 // dQ/dl = qc/(dx * sqrt(1 + dy/dx^2 + dz/dx^2))
189 //
190 // The following effects are included in the calculation:
191 // 1. calibration values for t0 and vdrift (using x coordinate to calculate slice)
192 // 2. cluster sharing (optional see AliTRDrecoParam::SetClusterSharing())
193 // 3. cluster size
194 //
195
196         Int_t nclusters[knSlices];
197         for(int i=0; i<knSlices; i++){ 
198                 fdEdx[i]     = 0.;
199                 nclusters[i] = 0;
200         }
201         Float_t clength = (/*.5 * */AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
202
203         AliTRDcluster *cluster = 0x0;
204         for(int ic=0; ic<fgTimeBins; ic++){
205                 if(!(cluster = fClusters[ic])) continue;
206                 Float_t x = cluster->GetX();
207                 
208                 // Filter clusters for dE/dx calculation
209                 
210                 // 1.consider calibration effects for slice determination
211                 Int_t slice; 
212                 if(cluster->IsInChamber()) slice = Int_t(TMath::Abs(fX0 - x) * nslices / clength);
213                 else slice = x < fX0 ? 0 : nslices-1;
214                 
215                 // 2. take sharing into account
216                 Float_t w = cluster->IsShared() ? .5 : 1.;
217                 
218                 // 3. take into account large clusters TODO
219                 //w *= c->GetNPads() > 3 ? .8 : 1.;
220                 
221                 //CHECK !!!
222                 fdEdx[slice]   += w * GetdQdl(ic); //fdQdl[ic];
223                 nclusters[slice]++;
224         } // End of loop over clusters
225
226         // calculate mean charge per slice
227         for(int is=0; is<nslices; is++) if(nclusters[is]) fdEdx[is] /= nclusters[is];
228 }
229
230 //____________________________________________________________________
231 Float_t AliTRDseedV1::GetdQdl(Int_t ic) const
232 {
233         return fClusters[ic] ? TMath::Abs(fClusters[ic]->GetQ()) /fdX / TMath::Sqrt(1. + fYfit[1]*fYfit[1] + fZfit[1]*fZfit[1]) : 0.;
234 }
235
236 //____________________________________________________________________
237 Double_t* AliTRDseedV1::GetProbability()
238 {       
239 // Fill probability array for tracklet from the DB.
240 //
241 // Parameters
242 //
243 // Output
244 //   returns pointer to the probability array and 0x0 if missing DB access 
245 //
246 // Detailed description
247
248         
249         // retrive calibration db
250   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
251   if (!calibration) {
252     AliError("No access to calibration data");
253     return 0x0;
254   }
255
256   // Retrieve the CDB container class with the parametric detector response
257   const AliTRDCalPID *pd = calibration->GetPIDObject(fRecoParam->GetPIDMethod());
258   if (!pd) {
259     AliError("No access to AliTRDCalPID object");
260     return 0x0;
261   }
262         
263         // calculate tracklet length TO DO
264   Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
265   /// TMath::Sqrt((1.0 - fSnp[iPlane]*fSnp[iPlane]) / (1.0 + fTgl[iPlane]*fTgl[iPlane]));
266   
267   //calculate dE/dx
268   CookdEdx(fRecoParam->GetNdEdxSlices());
269   
270   // Sets the a priori probabilities
271   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
272     fProb[ispec] = pd->GetProbability(ispec, fMom, &fdEdx[0], length, fPlane);  
273   }
274
275         return &fProb[0];
276 }
277
278 //____________________________________________________________________
279 Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
280 {
281   //
282   // Returns a quality measurement of the current seed
283   //
284
285         Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
286         return .5 * (18.0 - fN2)
287                 + 10.* TMath::Abs(fYfit[1] - fYref[1])
288                 + 5.* TMath::Abs(fYfit[0] - fYref[0] + zcorr)
289                 + 2. * TMath::Abs(fMeanz - fZref[0]) / fPadLength;
290 }
291
292 //____________________________________________________________________
293 void AliTRDseedV1::GetCovAt(Double_t /*x*/, Double_t *cov) const
294 {
295 // Computes covariance in the y-z plane at radial point x
296
297         const Float_t k0= .2; // to be checked in FindClusters
298         Double_t sy20   = k0*TMath::Tan(fYfit[1]); sy20 *= sy20;
299         
300         Double_t sy2    = fSigmaY2*fSigmaY2 + sy20;
301         Double_t sz2    = fPadLength/12.;
302
303         //printf("Yfit[1] %f sy20 %f SigmaY2 %f\n", fYfit[1], sy20, fSigmaY2);
304
305         cov[0] = sy2;
306         cov[1] = fTilt*(sy2-sz2);
307         cov[2] = sz2;
308 }
309
310
311 //____________________________________________________________________
312 void AliTRDseedV1::SetOwner(Bool_t own)
313 {
314         //AliInfo(Form("own [%s] fOwner[%s]", own?"YES":"NO", fOwner?"YES":"NO"));
315         
316         if(own){
317                 for(int ic=0; ic<knTimebins; ic++){
318                         if(!fClusters[ic]) continue;
319                         fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
320                 }
321                 fOwner = kTRUE;
322         } else {
323                 if(fOwner){
324                         for(int ic=0; ic<knTimebins; ic++){
325                                 if(!fClusters[ic]) continue;
326                                 delete fClusters[ic];
327                                 //fClusters[ic] = tracker->GetClusters(index) TODO
328                         }
329                 }
330                 fOwner = kFALSE;
331         }
332 }
333
334 //____________________________________________________________________
335 Bool_t  AliTRDseedV1::AttachClustersIter(AliTRDstackLayer *layer
336                                        , Float_t quality
337                                        , Bool_t kZcorr
338                                        , AliTRDcluster *c)
339 {
340   //
341   // Iterative process to register clusters to the seed.
342   // In iteration 0 we try only one pad-row and if quality not
343   // sufficient we try 2 pad-rows (about 5% of tracks cross 2 pad-rows)
344   //
345         
346         if(!fRecoParam){
347                 AliError("Seed can not be used without a valid RecoParam.");
348                 return kFALSE;
349         }
350
351         //AliInfo(Form("TimeBins = %d TimeBinsRange = %d", fgTimeBins, fTimeBinsRange));
352
353         Float_t  tquality;
354         Double_t kroady = fRecoParam->GetRoad1y();
355         Double_t kroadz = fPadLength * .5 + 1.;
356         
357         // initialize configuration parameters
358         Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
359         Int_t   niter = kZcorr ? 1 : 2;
360         
361         Double_t yexp, zexp;
362         Int_t ncl = 0;
363         // start seed update
364         for (Int_t iter = 0; iter < niter; iter++) {
365         //AliInfo(Form("iter = %i", iter));
366                 ncl = 0;
367                 for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
368                         // define searching configuration
369                         Double_t dxlayer = layer[iTime].GetX() - fX0;
370                         if(c){
371                                 zexp = c->GetZ();
372                                 //Try 2 pad-rows in second iteration
373                                 if (iter > 0) {
374                                         zexp = fZref[0] + fZref[1] * dxlayer - zcorr;
375                                         if (zexp > c->GetZ()) zexp = c->GetZ() + fPadLength*0.5;
376                                         if (zexp < c->GetZ()) zexp = c->GetZ() - fPadLength*0.5;
377                                 }
378                         } else zexp = fZref[0];
379                         yexp  = fYref[0] + fYref[1] * dxlayer - zcorr;
380                         
381                         // Get and register cluster
382                         Int_t    index = layer[iTime].SearchNearestCluster(yexp, zexp, kroady, kroadz);
383                         if (index < 0) continue;
384                         AliTRDcluster *cl = (AliTRDcluster*) layer[iTime].GetCluster(index);
385                         
386                         Int_t globalIndex = layer[iTime].GetGlobalIndex(index);
387                         fIndexes[iTime]  = globalIndex;
388                         fClusters[iTime] = cl;
389                         fY[iTime]        = cl->GetY();
390                         fZ[iTime]        = cl->GetZ();
391                         ncl++;
392                 }
393                 
394                 if(ncl){        
395                         // calculate length of the time bin (calibration aware)
396                         Int_t irp = 0; Float_t x[2]; Int_t tb[2];
397                         for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
398                                 if(!fClusters[iTime]) continue;
399                                 x[irp]  = fClusters[iTime]->GetX();
400                                 tb[irp] = iTime;
401                                 irp++;
402                                 if(irp==2) break;
403                         } 
404                         fdX = (x[1] - x[0]) / (tb[0] - tb[1]);
405         
406                         // update X0 from the clusters (calibration/alignment aware)
407                         for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
408                                 if(!layer[iTime].IsT0()) continue;
409                                 if(fClusters[iTime]){ 
410                                         fX0 = fClusters[iTime]->GetX();
411                                         break;
412                                 } else { // we have to infere the position of the anode wire from the other clusters
413                                         for (Int_t jTime = iTime+1; jTime < fgTimeBins; jTime++) {
414                                                 if(!fClusters[jTime]) continue;
415                                                 fX0 = fClusters[jTime]->GetX() + fdX * (jTime - iTime);
416                                         }
417                                         break;
418                                 }
419                         }       
420                         
421                         // update YZ reference point
422                         // TODO
423                         
424                         // update x reference positions (calibration/alignment aware)
425                         for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
426                                 if(!fClusters[iTime]) continue;
427                                 fX[iTime] = fClusters[iTime]->GetX() - fX0;
428                         } 
429                         
430                         AliTRDseed::Update();
431                 }
432                 
433                 if(IsOK()){
434                         tquality = GetQuality(kZcorr);
435                         if(tquality < quality) break;
436                         else quality = tquality;
437                 }
438                 kroadz *= 2.;
439         } // Loop: iter
440         if (!IsOK()) return kFALSE;
441
442         CookLabels();
443         UpdateUsed();
444         return kTRUE;   
445 }
446
447 //____________________________________________________________________
448 Bool_t  AliTRDseedV1::AttachClusters(AliTRDstackLayer *layer
449                                        ,Bool_t kZcorr)
450 {
451   //
452   // Projective algorithm to attach clusters to seeding tracklets
453   //
454   // Parameters
455   //
456   // Output
457   //
458   // Detailed description
459   // 1. Collapse x coordinate for the full detector plane
460   // 2. truncated mean on y (r-phi) direction
461   // 3. purge clusters
462   // 4. truncated mean on z direction
463   // 5. purge clusters
464   // 6. fit tracklet
465   //    
466
467         if(!fRecoParam){
468                 AliError("Seed can not be used without a valid RecoParam.");
469                 return kFALSE;
470         }
471
472         const Int_t kClusterCandidates = 2 * knTimebins;
473         
474         //define roads
475         Double_t kroady = fRecoParam->GetRoad1y();
476         Double_t kroadz = fPadLength * 1.5 + 1.;
477         // correction to y for the tilting angle
478         Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
479
480         // working variables
481         AliTRDcluster *clusters[kClusterCandidates];
482         Double_t cond[4], yexp[knTimebins], zexp[knTimebins],
483                 yres[kClusterCandidates], zres[kClusterCandidates];
484         Int_t ncl, *index = 0x0, tboundary[knTimebins];
485         
486         // Do cluster projection
487         Int_t nYclusters = 0; Bool_t kEXIT = kFALSE;
488         for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
489                 fX[iTime] = layer[iTime].GetX() - fX0;
490                 zexp[iTime] = fZref[0] + fZref[1] * fX[iTime];
491                 yexp[iTime] = fYref[0] + fYref[1] * fX[iTime] - zcorr;
492                 
493                 // build condition and process clusters
494                 cond[0] = yexp[iTime] - kroady; cond[1] = yexp[iTime] + kroady;
495                 cond[2] = zexp[iTime] - kroadz; cond[3] = zexp[iTime] + kroadz;
496                 layer[iTime].GetClusters(cond, index, ncl);
497                 for(Int_t ic = 0; ic<ncl; ic++){
498                         AliTRDcluster *c = layer[iTime].GetCluster(index[ic]);
499                         clusters[nYclusters] = c;
500                         yres[nYclusters++] = c->GetY() - yexp[iTime];
501                         if(nYclusters >= kClusterCandidates) {
502                                 AliWarning(Form("Cluster candidates reached limit %d. Some may be lost.", kClusterCandidates));
503                                 kEXIT = kTRUE;
504                                 break;
505                         }
506                 }
507                 tboundary[iTime] = nYclusters;
508                 if(kEXIT) break;
509         }
510         
511         // Evaluate truncated mean on the y direction
512         Double_t mean, sigma;
513         AliMathBase::EvaluateUni(nYclusters, yres, mean, sigma, Int_t(nYclusters*.8)-2);
514         //purge cluster candidates
515         Int_t nZclusters = 0;
516         for(Int_t ic = 0; ic<nYclusters; ic++){
517                 if(yres[ic] - mean > 4. * sigma){
518                         clusters[ic] = 0x0;
519                         continue;
520                 }
521                 zres[nZclusters++] = clusters[ic]->GetZ() - zexp[clusters[ic]->GetLocalTimeBin()];
522         }
523         
524         // Evaluate truncated mean on the z direction
525         AliMathBase::EvaluateUni(nZclusters, zres, mean, sigma, Int_t(nZclusters*.8)-2);
526         //purge cluster candidates
527         for(Int_t ic = 0; ic<nZclusters; ic++){
528                 if(zres[ic] - mean > 4. * sigma){
529                         clusters[ic] = 0x0;
530                         continue;
531                 }
532         }
533
534         
535         // Select only one cluster/TimeBin
536         Int_t lastCluster = 0;
537         fN2 = 0;
538         for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
539                 ncl = tboundary[iTime] - lastCluster;
540                 if(!ncl) continue;
541                 AliTRDcluster *c = 0x0;
542                 if(ncl == 1){
543                         c = clusters[lastCluster];
544                 } else if(ncl > 1){
545                         Float_t dold = 9999.; Int_t iptr = lastCluster;
546                         for(int ic=lastCluster; ic<tboundary[iTime]; ic++){
547                                 if(!clusters[ic]) continue;
548                                 Float_t y = yexp[iTime] - clusters[ic]->GetY();
549                                 Float_t z = zexp[iTime] - clusters[ic]->GetZ();
550                                 Float_t d = y * y + z * z;
551                                 if(d > dold) continue;
552                                 dold = d;
553                                 iptr = ic;
554                         }
555                         c = clusters[iptr];
556                 }
557                 //Int_t GlobalIndex = layer[iTime].GetGlobalIndex(index);
558                 //fIndexes[iTime]  = GlobalIndex;
559                 fClusters[iTime] = c;
560                 fY[iTime]        = c->GetY();
561                 fZ[iTime]        = c->GetZ();
562                 lastCluster      = tboundary[iTime];
563                 fN2++;
564         }
565         
566         // number of minimum numbers of clusters expected for the tracklet
567         Int_t kClmin = Int_t(fRecoParam->GetFindableClusters()*fgTimeBins);
568   if (fN2 < kClmin){
569                 AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin));
570     fN2 = 0;
571     return kFALSE;
572   }
573
574         // update used clusters
575         fNUsed = 0;
576         for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
577                 if(!fClusters[iTime]) continue;
578                 if((fClusters[iTime]->IsUsed())) fNUsed++;
579         }
580
581   if (fN2-fNUsed < kClmin){
582                 AliWarning(Form("Too many clusters already in use %d (from %d).", fNUsed, fN2));
583     fN2 = 0;
584     return kFALSE;
585   }
586         
587         return kTRUE;
588 }
589
590 //____________________________________________________________________
591 Bool_t AliTRDseedV1::Fit()
592 {
593   //
594   // Linear fit of the tracklet
595   //
596   // Parameters :
597   //
598   // Output :
599   //  True if successful
600   //
601   // Detailed description
602   // 2. Check if tracklet crosses pad row boundary
603   // 1. Calculate residuals in the y (r-phi) direction
604   // 3. Do a Least Square Fit to the data
605   //
606
607         //Float_t  sigmaexp  = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction
608   Float_t  ycrosscor = fPadLength * fTilt * 0.5; // Y correction for crossing
609   Float_t  anglecor = fTilt * fZref[1];  // Correction to the angle
610
611         // calculate residuals
612         Float_t yres[knTimebins]; // y (r-phi) residuals
613         Int_t zint[knTimebins],   // Histograming of the z coordinate
614               zout[2*knTimebins];//
615         
616         fN = 0;
617         for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
618     if (!fClusters[iTime]) continue;
619     if (!fClusters[iTime]->IsInChamber()) continue;
620     yres[iTime] = fY[iTime] - fYref[0] - (fYref[1] + anglecor) * fX[iTime] + fTilt * (fZ[iTime] - fZref[0]);
621                 zint[fN] = Int_t(fZ[iTime]);
622                 fN++;
623         }
624
625         // calculate pad row boundary crosses
626         Int_t kClmin = Int_t(fRecoParam->GetFindableClusters()*fgTimeBins);
627         Int_t nz = AliMathBase::Freq(fN, zint, zout, kFALSE);
628   fZProb   = zout[0];
629   if(nz <= 1) zout[3] = 0;
630   if(zout[1] + zout[3] < kClmin) {
631                 AliWarning(Form("Not enough clusters to fit the cross boundary tracklet %d [%d].", zout[1]+zout[3], kClmin));
632                 return kFALSE;
633         }
634   // Z distance bigger than pad - length
635   if (TMath::Abs(zout[0]-zout[2]) > fPadLength) zout[3]=0;
636  
637
638   Double_t sumw   = 0., 
639                 sumwx  = 0.,
640                 sumwx2 = 0.,
641                 sumwy  = 0.,
642                 sumwxy = 0.,
643                 sumwz  = 0.,
644                 sumwxz = 0.;
645         Int_t npads;
646   fMPads = 0;
647         fMeanz = 0.;
648         // we will use only the clusters which are in the detector range
649         for(int iTime=0; iTime<fgTimeBins; iTime++){
650     fUsable[iTime] = kFALSE;
651     if (!fClusters[iTime]) continue;
652                 npads = fClusters[iTime]->GetNPads();
653
654                 fUsable[iTime] = kTRUE;
655     fN2++;
656     fMPads += npads;
657     Float_t weight = 1.0;
658     if(npads > 5) weight = 0.2;
659     else if(npads > 4) weight = 0.5;
660     sumw   += weight; 
661     sumwx  += fX[iTime] * weight;
662     sumwx2 += fX[iTime] * fX[iTime] * weight;
663     sumwy  += weight * yres[iTime];
664     sumwxy += weight * yres[iTime] * fX[iTime];
665     sumwz  += weight * fZ[iTime];    
666     sumwxz += weight * fZ[iTime] * fX[iTime];
667         }
668   if (fN2 < kClmin){
669                 AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin));
670     fN2 = 0;
671     return kFALSE;
672   }
673   fMeanz = sumwz / sumw;
674         fNChange = 0;
675
676         // Tracklet on boundary
677   Float_t correction = 0;
678   if (fNChange > 0) {
679     if (fMeanz < fZProb) correction =  ycrosscor;
680     if (fMeanz > fZProb) correction = -ycrosscor;
681   }
682
683   Double_t det = sumw * sumwx2 - sumwx * sumwx;
684   fYfitR[0]    = (sumwx2 * sumwy  - sumwx * sumwxy) / det;
685   fYfitR[1]    = (sumw   * sumwxy - sumwx * sumwy)  / det;
686   
687   fSigmaY2 = 0;
688   for (Int_t i = 0; i < fgTimeBins+1; i++) {
689     if (!fUsable[i]) continue;
690     Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i];
691     fSigmaY2 += delta*delta;
692   }
693   fSigmaY2 = TMath::Sqrt(fSigmaY2 / Float_t(fN2-2));
694   
695   fZfitR[0]  = (sumwx2 * sumwz  - sumwx * sumwxz) / det;
696   fZfitR[1]  = (sumw   * sumwxz - sumwx * sumwz)  / det;
697   fZfit[0]   = (sumwx2 * sumwz  - sumwx * sumwxz) / det;
698   fZfit[1]   = (sumw   * sumwxz - sumwx * sumwz)  / det;
699   fYfitR[0] += fYref[0] + correction;
700   fYfitR[1] += fYref[1];
701   fYfit[0]   = fYfitR[0];
702   fYfit[1]   = fYfitR[1];
703
704         return kTRUE;
705 }
706
707 //_____________________________________________________________________________
708 Float_t AliTRDseedV1::FitRiemanTilt(AliTRDseedV1 *cseed, Bool_t terror)
709 {
710   //
711   // Fit the Rieman tilt
712   //
713
714   // Fitting with tilting pads - kz not fixed
715         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
716         Int_t nTimeBins = cal->GetNumberOfTimeBins();
717   TLinearFitter fitterT2(4,"hyp4");  
718   fitterT2.StoreData(kTRUE);
719   Float_t xref2 = (cseed[2].fX0 + cseed[3].fX0) * 0.5; // Reference x0 for z
720   
721   Int_t npointsT = 0;
722   fitterT2.ClearPoints();
723
724   for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
725 //              printf("\nLayer %d\n", iLayer);
726 //     cseed[iLayer].Print();
727                 if (!cseed[iLayer].IsOK()) continue;
728     Double_t tilt = cseed[iLayer].fTilt;
729
730     for (Int_t itime = 0; itime < nTimeBins+1; itime++) {
731 //                      printf("\ttime %d\n", itime);
732       if (!cseed[iLayer].fUsable[itime]) continue;
733       // x relative to the midle chamber
734       Double_t x = cseed[iLayer].fX[itime] + cseed[iLayer].fX0 - xref2;  
735       Double_t y = cseed[iLayer].fY[itime];
736       Double_t z = cseed[iLayer].fZ[itime];
737
738       //
739       // Tilted rieman
740       //
741       Double_t uvt[6];
742       Double_t x2 = cseed[iLayer].fX[itime] + cseed[iLayer].fX0;      // Global x
743       Double_t t  = 1.0 / (x2*x2 + y*y);
744       uvt[1]  = t;
745       uvt[0]  = 2.0 * x2   * uvt[1];
746       uvt[2]  = 2.0 * tilt * uvt[1];
747       uvt[3]  = 2.0 * tilt *uvt[1] * x;       
748       uvt[4]  = 2.0 * (y + tilt * z) * uvt[1];
749       
750       Double_t error = 2.0 * uvt[1];
751       error *= terror ? cseed[iLayer].fSigmaY : .2;
752
753 //                      printf("\tadd point :\n");
754 //                      for(int i=0; i<5; i++) printf("%f ", uvt[i]);
755 //                      printf("\n");
756       fitterT2.AddPoint(uvt,uvt[4],error);
757       npointsT++;
758
759     }
760
761   }
762   fitterT2.Eval();
763   Double_t rpolz0 = fitterT2.GetParameter(3);
764   Double_t rpolz1 = fitterT2.GetParameter(4);       
765
766   //
767   // Linear fitter  - not possible to make boundaries
768   // non accept non possible z and dzdx combination
769   //        
770   Bool_t acceptablez = kTRUE;
771   for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
772     if (cseed[iLayer].IsOK()) {
773       Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].fX0 - xref2);
774       if (TMath::Abs(cseed[iLayer].fZProb - zT2) > cseed[iLayer].fPadLength * 0.5 + 1.0) {
775         acceptablez = kFALSE;
776       }
777     }
778   }
779   if (!acceptablez) {
780     Double_t zmf  = cseed[2].fZref[0] + cseed[2].fZref[1] * (xref2 - cseed[2].fX0);
781     Double_t dzmf = (cseed[2].fZref[1] + cseed[3].fZref[1]) * 0.5;
782     fitterT2.FixParameter(3,zmf);
783     fitterT2.FixParameter(4,dzmf);
784     fitterT2.Eval();
785     fitterT2.ReleaseParameter(3);
786     fitterT2.ReleaseParameter(4);
787     rpolz0 = fitterT2.GetParameter(3);
788     rpolz1 = fitterT2.GetParameter(4);
789   }
790   
791   Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT);  
792   Double_t params[3];
793   params[0] =  fitterT2.GetParameter(0);
794   params[1] =  fitterT2.GetParameter(1);
795   params[2] =  fitterT2.GetParameter(2);            
796   Double_t curvature =  1.0 + params[1] * params[1] - params[2] * params[0];
797
798   for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
799
800     Double_t  x  = cseed[iLayer].fX0;
801     Double_t  y  = 0;
802     Double_t  dy = 0;
803     Double_t  z  = 0;
804     Double_t  dz = 0;
805
806     // y
807     Double_t res2 = (x * params[0] + params[1]);
808     res2 *= res2;
809     res2  = 1.0 - params[2]*params[0] + params[1]*params[1] - res2;
810     if (res2 >= 0) {
811       res2 = TMath::Sqrt(res2);
812       y    = (1.0 - res2) / params[0];
813     }
814
815     //dy
816     Double_t x0 = -params[1] / params[0];
817     if (-params[2]*params[0] + params[1]*params[1] + 1 > 0) {
818       Double_t rm1 = params[0] / TMath::Sqrt(-params[2]*params[0] + params[1]*params[1] + 1); 
819       if (1.0/(rm1*rm1) - (x-x0) * (x-x0) > 0.0) {
820                                 Double_t res = (x - x0) / TMath::Sqrt(1.0 / (rm1*rm1) - (x-x0)*(x-x0));
821                                 if (params[0] < 0) res *= -1.0;
822                                 dy = res;
823       }
824     }
825     z  = rpolz0 + rpolz1 * (x - xref2);
826     dz = rpolz1;
827     cseed[iLayer].fYref[0] = y;
828     cseed[iLayer].fYref[1] = dy;
829     cseed[iLayer].fZref[0] = z;
830     cseed[iLayer].fZref[1] = dz;
831     cseed[iLayer].fC       = curvature;
832     
833   }
834
835   return chi2TR;
836
837 }
838
839 //___________________________________________________________________
840 void AliTRDseedV1::Print()
841 {
842   //
843   // Printing the seedstatus
844   //
845
846         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
847         Int_t nTimeBins = cal->GetNumberOfTimeBins();
848         
849         printf("Seed status :\n");
850         printf("  fTilt      = %f\n", fTilt);
851         printf("  fPadLength = %f\n", fPadLength);
852         printf("  fX0        = %f\n", fX0);
853         for(int ic=0; ic<nTimeBins; ic++) {
854           const Char_t *isUsable = fUsable[ic]?"Yes":"No";
855           printf("  %d X[%f] Y[%f] Z[%f] Indexes[%d] clusters[%p] usable[%s]\n"
856                 , ic
857                 , fX[ic]
858                 , fY[ic]
859                 , fZ[ic]
860                 , fIndexes[ic]
861                 , ((void*) fClusters[ic])
862                 , isUsable);
863         }
864
865         printf("  fYref[0] =%f fYref[1] =%f\n", fYref[0], fYref[1]);
866         printf("  fZref[0] =%f fZref[1] =%f\n", fZref[0], fZref[1]);
867         printf("  fYfit[0] =%f fYfit[1] =%f\n", fYfit[0], fYfit[1]);
868         printf("  fYfitR[0]=%f fYfitR[1]=%f\n", fYfitR[0], fYfitR[1]);
869         printf("  fZfit[0] =%f fZfit[1] =%f\n", fZfit[0], fZfit[1]);
870         printf("  fZfitR[0]=%f fZfitR[1]=%f\n", fZfitR[0], fZfitR[1]);
871         printf("  fSigmaY =%f\n", fSigmaY);
872         printf("  fSigmaY2=%f\n", fSigmaY2);            
873         printf("  fMeanz  =%f\n", fMeanz);
874         printf("  fZProb  =%f\n", fZProb);
875         printf("  fLabels[0]=%d fLabels[1]=%d\n", fLabels[0], fLabels[1]);
876         printf("  fN      =%d\n", fN);
877         printf("  fN2     =%d (>8 isOK)\n",fN2);
878         printf("  fNUsed  =%d\n", fNUsed);
879         printf("  fFreq   =%d\n", fFreq);
880         printf("  fNChange=%d\n",  fNChange);
881         printf("  fMPads  =%f\n", fMPads);
882         
883         printf("  fC      =%f\n", fC);        
884         printf("  fCC     =%f\n",fCC);      
885         printf("  fChi2   =%f\n", fChi2);  
886         printf("  fChi2Z  =%f\n", fChi2Z);
887
888 }