43792bc2278996ce31e336f708ab6c64eb71a164
[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 #include "TClonesArray.h" // tmp
31 #include <TTreeStream.h>
32
33 #include "AliLog.h"
34 #include "AliMathBase.h"
35
36 #include "AliTRDcluster.h"
37 #include "AliTRDseedV1.h"
38 #include "AliTRDtrackV1.h"
39 #include "AliTRDcalibDB.h"
40 #include "AliTRDchamberTimeBin.h"
41 #include "AliTRDtrackingChamber.h"
42 #include "AliTRDtrackerV1.h"
43 #include "AliTRDReconstructor.h"
44 #include "AliTRDrecoParam.h"
45 #include "Cal/AliTRDCalPID.h"
46
47 ClassImp(AliTRDseedV1)
48
49 //____________________________________________________________________
50 AliTRDseedV1::AliTRDseedV1(Int_t det) 
51   :AliTRDseed()
52   ,fReconstructor(0x0)
53   ,fClusterIter(0x0)
54   ,fClusterIdx(0)
55   ,fDet(det)
56   ,fMom(0.)
57   ,fSnp(0.)
58   ,fTgl(0.)
59   ,fdX(0.)
60 {
61   //
62   // Constructor
63   //
64   //printf("AliTRDseedV1::AliTRDseedV1()\n");
65
66   for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = 0.;
67   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec]  = -1.;
68 }
69
70 //____________________________________________________________________
71 AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
72   :AliTRDseed((AliTRDseed&)ref)
73   ,fReconstructor(ref.fReconstructor)
74   ,fClusterIter(0x0)
75   ,fClusterIdx(0)
76   ,fDet(ref.fDet)
77   ,fMom(ref.fMom)
78   ,fSnp(ref.fSnp)
79   ,fTgl(ref.fTgl)
80   ,fdX(ref.fdX)
81 {
82   //
83   // Copy Constructor performing a deep copy
84   //
85
86   //printf("AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &)\n");
87   SetBit(kOwner, kFALSE);
88   for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = ref.fdEdx[islice];
89   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = ref.fProb[ispec];
90 }
91
92
93 //____________________________________________________________________
94 AliTRDseedV1& AliTRDseedV1::operator=(const AliTRDseedV1 &ref)
95 {
96   //
97   // Assignment Operator using the copy function
98   //
99
100   if(this != &ref){
101     ref.Copy(*this);
102   }
103   return *this;
104
105 }
106
107 //____________________________________________________________________
108 AliTRDseedV1::~AliTRDseedV1()
109 {
110   //
111   // Destructor. The RecoParam object belongs to the underlying tracker.
112   //
113
114   //printf("I-AliTRDseedV1::~AliTRDseedV1() : Owner[%s]\n", IsOwner()?"YES":"NO");
115
116   if(IsOwner()) 
117     for(int itb=0; itb<knTimebins; itb++){
118       if(!fClusters[itb]) continue; 
119       //AliInfo(Form("deleting c %p @ %d", fClusters[itb], itb));
120       delete fClusters[itb];
121       fClusters[itb] = 0x0;
122     }
123 }
124
125 //____________________________________________________________________
126 void AliTRDseedV1::Copy(TObject &ref) const
127 {
128   //
129   // Copy function
130   //
131
132   //AliInfo("");
133   AliTRDseedV1 &target = (AliTRDseedV1 &)ref; 
134
135   target.fClusterIter   = 0x0;
136   target.fClusterIdx    = 0;
137   target.fDet           = fDet;
138   target.fMom           = fMom;
139   target.fSnp           = fSnp;
140   target.fTgl           = fTgl;
141   target.fdX            = fdX;
142   target.fReconstructor = fReconstructor;
143   
144   for(int islice=0; islice < knSlices; islice++) target.fdEdx[islice] = fdEdx[islice];
145   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) target.fProb[ispec] = fProb[ispec];
146   
147   AliTRDseed::Copy(target);
148 }
149
150
151 //____________________________________________________________
152 Bool_t AliTRDseedV1::Init(AliTRDtrackV1 *track)
153 {
154 // Initialize this tracklet using the track information
155 //
156 // Parameters:
157 //   track - the TRD track used to initialize the tracklet
158 // 
159 // Detailed description
160 // The function sets the starting point and direction of the
161 // tracklet according to the information from the TRD track.
162 // 
163 // Caution
164 // The TRD track has to be propagated to the beginning of the
165 // chamber where the tracklet will be constructed
166 //
167
168   Double_t y, z; 
169   if(!track->GetProlongation(fX0, y, z)) return kFALSE;
170   fYref[0] = y;
171   fYref[1] = track->GetSnp()/(1. - track->GetSnp()*track->GetSnp());
172   fZref[0] = z;
173   fZref[1] = track->GetTgl();
174
175   //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());
176   return kTRUE;
177 }
178
179
180 //____________________________________________________________________
181 void AliTRDseedV1::CookdEdx(Int_t nslices)
182 {
183 // Calculates average dE/dx for all slices and store them in the internal array fdEdx. 
184 //
185 // Parameters:
186 //  nslices : number of slices for which dE/dx should be calculated
187 // Output:
188 //  store results in the internal array fdEdx. This can be accessed with the method
189 //  AliTRDseedV1::GetdEdx()
190 //
191 // Detailed description
192 // Calculates average dE/dx for all slices. Depending on the PID methode 
193 // the number of slices can be 3 (LQ) or 8(NN). 
194 // The calculation of dQ/dl are done using the tracklet fit results (see AliTRDseedV1::GetdQdl(Int_t)) i.e.
195 //
196 // dQ/dl = qc/(dx * sqrt(1 + dy/dx^2 + dz/dx^2))
197 //
198 // The following effects are included in the calculation:
199 // 1. calibration values for t0 and vdrift (using x coordinate to calculate slice)
200 // 2. cluster sharing (optional see AliTRDrecoParam::SetClusterSharing())
201 // 3. cluster size
202 //
203
204   Int_t nclusters[knSlices];
205   for(int i=0; i<knSlices; i++){ 
206     fdEdx[i]     = 0.;
207     nclusters[i] = 0;
208   }
209   Float_t clength = (/*.5 * */AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
210
211   AliTRDcluster *cluster = 0x0;
212   for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++){
213     if(!(cluster = fClusters[ic])) continue;
214     Float_t x = cluster->GetX();
215     
216     // Filter clusters for dE/dx calculation
217     
218     // 1.consider calibration effects for slice determination
219     Int_t slice; 
220     if(cluster->IsInChamber()) slice = Int_t(TMath::Abs(fX0 - x) * nslices / clength);
221     else slice = x < fX0 ? 0 : nslices-1;
222     
223     // 2. take sharing into account
224     Float_t w = cluster->IsShared() ? .5 : 1.;
225     
226     // 3. take into account large clusters TODO
227     //w *= c->GetNPads() > 3 ? .8 : 1.;
228     
229     //CHECK !!!
230     fdEdx[slice]   += w * GetdQdl(ic); //fdQdl[ic];
231     nclusters[slice]++;
232   } // End of loop over clusters
233
234   //if(fReconstructor->GetPIDMethod() == AliTRDReconstructor::kLQPID){
235   if(nslices == AliTRDReconstructor::kLQslices){
236   // calculate mean charge per slice (only LQ PID)
237     for(int is=0; is<nslices; is++){ 
238       if(nclusters[is]) fdEdx[is] /= nclusters[is];
239     }
240   }
241 }
242
243
244 //____________________________________________________________________
245 Float_t AliTRDseedV1::GetdQdl(Int_t ic) const
246 {
247   return fClusters[ic] ? TMath::Abs(fClusters[ic]->GetQ()) /fdX / TMath::Sqrt(1. + fYfit[1]*fYfit[1] + fZref[1]*fZref[1]) : 0.;
248 }
249
250 //____________________________________________________________________
251 Double_t* AliTRDseedV1::GetProbability()
252 {       
253 // Fill probability array for tracklet from the DB.
254 //
255 // Parameters
256 //
257 // Output
258 //   returns pointer to the probability array and 0x0 if missing DB access 
259 //
260 // Detailed description
261
262   
263   // retrive calibration db
264   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
265   if (!calibration) {
266     AliError("No access to calibration data");
267     return 0x0;
268   }
269
270   if (!fReconstructor) {
271     AliError("Reconstructor not set.");
272     return 0x0;
273   }
274
275   // Retrieve the CDB container class with the parametric detector response
276   const AliTRDCalPID *pd = calibration->GetPIDObject(fReconstructor->GetPIDMethod());
277   if (!pd) {
278     AliError("No access to AliTRDCalPID object");
279     return 0x0;
280   }
281   //AliInfo(Form("Method[%d] : %s", fReconstructor->GetRecoParam() ->GetPIDMethod(), pd->IsA()->GetName()));
282
283   // calculate tracklet length TO DO
284   Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
285   /// TMath::Sqrt((1.0 - fSnp[iPlane]*fSnp[iPlane]) / (1.0 + fTgl[iPlane]*fTgl[iPlane]));
286   
287   //calculate dE/dx
288   CookdEdx(fReconstructor->GetNdEdxSlices());
289   
290   // Sets the a priori probabilities
291   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
292     fProb[ispec] = pd->GetProbability(ispec, fMom, &fdEdx[0], length, GetPlane());      
293   }
294
295   return &fProb[0];
296 }
297
298 //____________________________________________________________________
299 Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
300 {
301   //
302   // Returns a quality measurement of the current seed
303   //
304
305   Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
306   return 
307       .5 * TMath::Abs(18.0 - fN2)
308     + 10.* TMath::Abs(fYfit[1] - fYref[1])
309     + 5. * TMath::Abs(fYfit[0] - fYref[0] + zcorr)
310     + 2. * TMath::Abs(fMeanz - fZref[0]) / fPadLength;
311 }
312
313 //____________________________________________________________________
314 void AliTRDseedV1::GetCovAt(Double_t /*x*/, Double_t *cov) const
315 {
316 // Computes covariance in the y-z plane at radial point x
317
318   Int_t ic = 0; while (!fClusters[ic]) ic++; 
319   AliTRDcalibDB *fCalib = AliTRDcalibDB::Instance();
320   Double_t exB         = fCalib->GetOmegaTau(fCalib->GetVdriftAverage(fClusters[ic]->GetDetector()), -AliTracker::GetBz()*0.1);
321
322   Double_t sy2    = fSigmaY2*fSigmaY2 + .2*(fYfit[1]-exB)*(fYfit[1]-exB);
323   Double_t sz2    = fPadLength/12.;
324
325
326   //printf("Yfit[1] %f sy20 %f SigmaY2 %f\n", fYfit[1], sy20, fSigmaY2);
327
328   cov[0] = sy2;
329   cov[1] = fTilt*(sy2-sz2);
330   cov[2] = sz2;
331
332   // insert systematic uncertainties calibration and misalignment
333   Double_t sys[15];
334   fReconstructor->GetRecoParam()->GetSysCovMatrix(sys);
335   cov[0] += (sys[0]*sys[0]);
336   cov[2] += (sys[1]*sys[1]);
337 }
338
339
340 //____________________________________________________________________
341 void AliTRDseedV1::SetOwner()
342 {
343   //AliInfo(Form("own [%s] fOwner[%s]", own?"YES":"NO", fOwner?"YES":"NO"));
344   
345   if(TestBit(kOwner)) return;
346   for(int ic=0; ic<knTimebins; ic++){
347     if(!fClusters[ic]) continue;
348     fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
349   }
350   SetBit(kOwner);
351 }
352
353 //____________________________________________________________________
354 Bool_t  AliTRDseedV1::AttachClustersIter(AliTRDtrackingChamber *chamber, Float_t quality, Bool_t kZcorr, AliTRDcluster *c)
355 {
356   //
357   // Iterative process to register clusters to the seed.
358   // In iteration 0 we try only one pad-row and if quality not
359   // sufficient we try 2 pad-rows (about 5% of tracks cross 2 pad-rows)
360   //
361   // debug level 7
362   //
363   
364   if(!fReconstructor->GetRecoParam() ){
365     AliError("Seed can not be used without a valid RecoParam.");
366     return kFALSE;
367   }
368
369   AliTRDchamberTimeBin *layer = 0x0;
370   if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7){
371     AliTRDtrackingChamber *ch = new AliTRDtrackingChamber(*chamber); 
372     (*AliTRDtrackerV1::DebugStreamer()) << "AttachClustersIter"
373       << "chamber.="   << ch
374       << "tracklet.="  << this
375       << "\n";  
376   }
377
378   Float_t  tquality;
379   Double_t kroady = fReconstructor->GetRecoParam() ->GetRoad1y();
380   Double_t kroadz = fPadLength * .5 + 1.;
381   
382   // initialize configuration parameters
383   Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
384   Int_t   niter = kZcorr ? 1 : 2;
385   
386   Double_t yexp, zexp;
387   Int_t ncl = 0;
388   // start seed update
389   for (Int_t iter = 0; iter < niter; iter++) {
390     ncl = 0;
391     for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
392       if(!(layer = chamber->GetTB(iTime))) continue;
393       if(!Int_t(*layer)) continue;
394       
395       // define searching configuration
396       Double_t dxlayer = layer->GetX() - fX0;
397       if(c){
398         zexp = c->GetZ();
399         //Try 2 pad-rows in second iteration
400         if (iter > 0) {
401           zexp = fZref[0] + fZref[1] * dxlayer - zcorr;
402           if (zexp > c->GetZ()) zexp = c->GetZ() + fPadLength*0.5;
403           if (zexp < c->GetZ()) zexp = c->GetZ() - fPadLength*0.5;
404         }
405       } else zexp = fZref[0] + (kZcorr ? fZref[1] * dxlayer : 0.);
406       yexp  = fYref[0] + fYref[1] * dxlayer - zcorr;
407       
408       // Get and register cluster
409       Int_t    index = layer->SearchNearestCluster(yexp, zexp, kroady, kroadz);
410       if (index < 0) continue;
411       AliTRDcluster *cl = (*layer)[index];
412       
413       fIndexes[iTime]  = layer->GetGlobalIndex(index);
414       fClusters[iTime] = cl;
415       fY[iTime]        = cl->GetY();
416       fZ[iTime]        = cl->GetZ();
417       ncl++;
418     }
419     if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d ncl [%d] = %d", iter, fDet, ncl));
420     
421     if(ncl>1){  
422       // calculate length of the time bin (calibration aware)
423       Int_t irp = 0; Float_t x[2]; Int_t tb[2];
424       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
425         if(!fClusters[iTime]) continue;
426         x[irp]  = fClusters[iTime]->GetX();
427         tb[irp] = iTime;
428         irp++;
429         if(irp==2) break;
430       } 
431       fdX = (x[1] - x[0]) / (tb[0] - tb[1]);
432   
433       // update X0 from the clusters (calibration/alignment aware)
434       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
435         if(!(layer = chamber->GetTB(iTime))) continue;
436         if(!layer->IsT0()) continue;
437         if(fClusters[iTime]){ 
438           fX0 = fClusters[iTime]->GetX();
439           break;
440         } else { // we have to infere the position of the anode wire from the other clusters
441           for (Int_t jTime = iTime+1; jTime < AliTRDtrackerV1::GetNTimeBins(); jTime++) {
442             if(!fClusters[jTime]) continue;
443             fX0 = fClusters[jTime]->GetX() + fdX * (jTime - iTime);
444           }
445           break;
446         }
447       } 
448       
449       // update YZ reference point
450       // TODO
451       
452       // update x reference positions (calibration/alignment aware)
453       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
454         if(!fClusters[iTime]) continue;
455         fX[iTime] = fClusters[iTime]->GetX() - fX0;
456       } 
457       
458       AliTRDseed::Update();
459     }
460     if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d nclFit [%d] = %d", iter, fDet, fN2));
461     
462     if(IsOK()){
463       tquality = GetQuality(kZcorr);
464       if(tquality < quality) break;
465       else quality = tquality;
466     }
467     kroadz *= 2.;
468   } // Loop: iter
469   if (!IsOK()) return kFALSE;
470
471   if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=1) CookLabels();
472   UpdateUsed();
473   return kTRUE; 
474 }
475
476 //____________________________________________________________________
477 Bool_t  AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *chamber
478                                       ,Bool_t kZcorr)
479 {
480   //
481   // Projective algorithm to attach clusters to seeding tracklets
482   //
483   // Parameters
484   //
485   // Output
486   //
487   // Detailed description
488   // 1. Collapse x coordinate for the full detector plane
489   // 2. truncated mean on y (r-phi) direction
490   // 3. purge clusters
491   // 4. truncated mean on z direction
492   // 5. purge clusters
493   // 6. fit tracklet
494   //    
495
496   if(!fReconstructor->GetRecoParam() ){
497     AliError("Seed can not be used without a valid RecoParam.");
498     return kFALSE;
499   }
500
501   const Int_t kClusterCandidates = 2 * knTimebins;
502   
503   //define roads
504   Double_t kroady = fReconstructor->GetRecoParam() ->GetRoad1y();
505   Double_t kroadz = fPadLength * 1.5 + 1.;
506   // correction to y for the tilting angle
507   Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
508
509   // working variables
510   AliTRDcluster *clusters[kClusterCandidates];
511   Double_t cond[4], yexp[knTimebins], zexp[knTimebins],
512     yres[kClusterCandidates], zres[kClusterCandidates];
513   Int_t ncl, *index = 0x0, tboundary[knTimebins];
514   
515   // Do cluster projection
516   AliTRDchamberTimeBin *layer = 0x0;
517   Int_t nYclusters = 0; Bool_t kEXIT = kFALSE;
518   for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
519     if(!(layer = chamber->GetTB(iTime))) continue;
520     if(!Int_t(*layer)) continue;
521     
522     fX[iTime] = layer->GetX() - fX0;
523     zexp[iTime] = fZref[0] + fZref[1] * fX[iTime];
524     yexp[iTime] = fYref[0] + fYref[1] * fX[iTime] - zcorr;
525     
526     // build condition and process clusters
527     cond[0] = yexp[iTime] - kroady; cond[1] = yexp[iTime] + kroady;
528     cond[2] = zexp[iTime] - kroadz; cond[3] = zexp[iTime] + kroadz;
529     layer->GetClusters(cond, index, ncl);
530     for(Int_t ic = 0; ic<ncl; ic++){
531       AliTRDcluster *c = layer->GetCluster(index[ic]);
532       clusters[nYclusters] = c;
533       yres[nYclusters++] = c->GetY() - yexp[iTime];
534       if(nYclusters >= kClusterCandidates) {
535         AliWarning(Form("Cluster candidates reached limit %d. Some may be lost.", kClusterCandidates));
536         kEXIT = kTRUE;
537         break;
538       }
539     }
540     tboundary[iTime] = nYclusters;
541     if(kEXIT) break;
542   }
543   
544   // Evaluate truncated mean on the y direction
545   Double_t mean, sigma;
546   AliMathBase::EvaluateUni(nYclusters, yres, mean, sigma, Int_t(nYclusters*.8)-2);
547   // purge cluster candidates
548   Int_t nZclusters = 0;
549   for(Int_t ic = 0; ic<nYclusters; ic++){
550     if(yres[ic] - mean > 4. * sigma){
551       clusters[ic] = 0x0;
552       continue;
553     }
554     zres[nZclusters++] = clusters[ic]->GetZ() - zexp[clusters[ic]->GetLocalTimeBin()];
555   }
556   
557   // Evaluate truncated mean on the z direction
558   AliMathBase::EvaluateUni(nZclusters, zres, mean, sigma, Int_t(nZclusters*.8)-2);
559   // purge cluster candidates
560   for(Int_t ic = 0; ic<nZclusters; ic++){
561     if(zres[ic] - mean > 4. * sigma){
562       clusters[ic] = 0x0;
563       continue;
564     }
565   }
566
567   
568   // Select only one cluster/TimeBin
569   Int_t lastCluster = 0;
570   fN2 = 0;
571   for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
572     ncl = tboundary[iTime] - lastCluster;
573     if(!ncl) continue;
574     Int_t iptr = lastCluster;
575     if(ncl > 1){
576       Float_t dold = 9999.;
577       for(int ic=lastCluster; ic<tboundary[iTime]; ic++){
578         if(!clusters[ic]) continue;
579         Float_t y = yexp[iTime] - clusters[ic]->GetY();
580         Float_t z = zexp[iTime] - clusters[ic]->GetZ();
581         Float_t d = y * y + z * z;
582         if(d > dold) continue;
583         dold = d;
584         iptr = ic;
585       }
586     }
587     fIndexes[iTime]  = chamber->GetTB(iTime)->GetGlobalIndex(iptr);
588     fClusters[iTime] = clusters[iptr];
589     fY[iTime]        = clusters[iptr]->GetY();
590     fZ[iTime]        = clusters[iptr]->GetZ();
591     lastCluster      = tboundary[iTime];
592     fN2++;
593   }
594   
595   // number of minimum numbers of clusters expected for the tracklet
596   Int_t kClmin = Int_t(fReconstructor->GetRecoParam() ->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
597   if (fN2 < kClmin){
598     AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin));
599     fN2 = 0;
600     return kFALSE;
601   }
602
603   // update used clusters
604   fNUsed = 0;
605   for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
606     if(!fClusters[iTime]) continue;
607     if((fClusters[iTime]->IsUsed())) fNUsed++;
608   }
609
610   if (fN2-fNUsed < kClmin){
611     AliWarning(Form("Too many clusters already in use %d (from %d).", fNUsed, fN2));
612     fN2 = 0;
613     return kFALSE;
614   }
615   
616   return kTRUE;
617 }
618
619 //____________________________________________________________________
620 Bool_t AliTRDseedV1::Fit(Bool_t tilt)
621 {
622   //
623   // Linear fit of the tracklet
624   //
625   // Parameters :
626   //
627   // Output :
628   //  True if successful
629   //
630   // Detailed description
631   // 2. Check if tracklet crosses pad row boundary
632   // 1. Calculate residuals in the y (r-phi) direction
633   // 3. Do a Least Square Fit to the data
634   //
635
636   const Int_t kClmin = 8;
637   const Float_t q0 = 100.;
638   const Float_t clSigma0 = 2.E-2;    //[cm]
639   const Float_t clSlopeQ = -1.19E-2; //[1/cm]
640
641   // get track direction
642   Double_t y0   = fYref[0];
643   Double_t dydx = fYref[1]; 
644   Double_t z0   = fZref[0];
645   Double_t dzdx = fZref[1];
646   Double_t yt, zt;
647
648   const Int_t kNtb = AliTRDtrackerV1::GetNTimeBins();
649   AliTRDtrackerV1::AliTRDLeastSquare fitterY, fitterZ;
650
651   // convertion factor from square to gauss distribution for sigma
652   Double_t convert = 1./TMath::Sqrt(12.);
653   
654   // book cluster information
655   Double_t xc[knTimebins], yc[knTimebins], zc[knTimebins], sy[knTimebins], sz[knTimebins];
656   Int_t zRow[knTimebins];
657
658
659   fN = 0;
660   AliTRDcluster *c=0x0, **jc = &fClusters[0];
661   for (Int_t ic=0; ic<kNtb; ic++, ++jc) {
662     zRow[ic] = -1;
663     xc[ic]  = -1.;
664     yc[ic]  = 999.;
665     zc[ic]  = 999.;
666     sy[ic]  = 0.;
667     sz[ic]  = 0.;
668     if(!(c = (*jc))) continue;
669     if(!c->IsInChamber()) continue;
670     Float_t w = 1.;
671     if(c->GetNPads()>4) w = .5;
672     if(c->GetNPads()>5) w = .2;
673     zRow[fN] = c->GetPadRow();
674     xc[fN]   = fX0 - c->GetX();
675     yc[fN]   = c->GetY();
676     zc[fN]   = c->GetZ();
677
678     // extrapolated y value for the track
679     yt = y0 - xc[fN]*dydx; 
680     // extrapolated z value for the track
681     zt = z0 - xc[fN]*dzdx; 
682     // tilt correction
683     if(tilt) yc[fN] -= fTilt*(zc[fN] - zt); 
684
685     // elaborate cluster error
686     Float_t qr = c->GetQ() - q0;
687     sy[fN]   = qr < 0. ? clSigma0*TMath::Exp(clSlopeQ*qr) : clSigma0;
688
689     fitterY.AddPoint(&xc[fN], yc[fN]-yt, sy[fN]);
690
691     sz[fN]   = fPadLength*convert;
692     fitterZ.AddPoint(&xc[fN], zc[fN], sz[fN]);
693     fN++;
694   }
695   // to few clusters
696   if (fN < kClmin) return kFALSE; 
697
698   // fit XY plane
699   fitterY.Eval();
700   fYfit[0] = y0+fitterY.GetFunctionParameter(0);
701   fYfit[1] = dydx-fitterY.GetFunctionParameter(1);
702
703   // check par row crossing
704   Int_t zN[2*AliTRDseed::knTimebins];
705   Int_t nz = AliTRDtrackerV1::Freq(fN, zRow, zN, kFALSE);
706   // more than one pad row crossing
707   if(nz>2) return kFALSE; 
708
709
710   // determine z offset of the fit
711   Float_t zslope = 0.;
712   Int_t nchanges = 0, nCross = 0;
713   if(nz==2){ // tracklet is crossing pad row
714     // Find the break time allowing one chage on pad-rows
715     // with maximal number of accepted clusters
716     Int_t padRef = zRow[0];
717     for (Int_t ic=1; ic<fN; ic++) {
718       if(zRow[ic] == padRef) continue;
719       
720       // debug
721       if(zRow[ic-1] == zRow[ic]){
722         printf("ERROR in pad row change!!!\n");
723       }
724     
725       // evaluate parameters of the crossing point
726       Float_t sx = (xc[ic-1] - xc[ic])*convert;
727       fCross[0] = .5 * (xc[ic-1] + xc[ic]);
728       fCross[2] = .5 * (zc[ic-1] + zc[ic]);
729       fCross[3] = TMath::Max(dzdx * sx, .01);
730       zslope    = zc[ic-1] > zc[ic] ? 1. : -1.;
731       padRef    = zRow[ic];
732       nCross    = ic;
733       nchanges++;
734     }
735   }
736
737   // condition on nCross and reset nchanges TODO
738
739   if(nchanges==1){
740     if(dzdx * zslope < 0.){
741       AliInfo("tracklet direction does not correspond to the track direction. TODO.");
742     }
743     SetBit(kRowCross, kTRUE); // mark pad row crossing
744     fitterZ.AddPoint(&fCross[0], fCross[2], fCross[3]);
745     fitterZ.Eval();
746     //zc[nc] = fitterZ.GetFunctionParameter(0); 
747     fCross[1] = fYfit[0] - fCross[0] * fYfit[1];
748     fCross[0] = fX0 - fCross[0];
749   } else if(nchanges > 1){ // debug
750     AliError("N pad row crossing > 1.");
751     return kFALSE;
752   }
753
754   UpdateUsed();
755
756   return kTRUE;
757 }
758
759
760 //___________________________________________________________________
761 void AliTRDseedV1::Print(Option_t*) const
762 {
763   //
764   // Printing the seedstatus
765   //
766
767   AliInfo(Form("Tracklet X0[%7.2f] Det[%d]", fX0, fDet));
768   printf("  Tilt[%+6.2f] PadLength[%5.2f]\n", fTilt, fPadLength);
769   AliTRDcluster* const* jc = &fClusters[0];
770   for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++, jc++) {
771     if(!(*jc)) continue;
772     printf("  %2d X[%7.2f] Y[%7.2f] Z[%7.2f] Idx[%d] c[%p] usable[%s]\n", 
773       ic, (*jc)->GetX(), (*jc)->GetY(), (*jc)->GetZ(), 
774       fIndexes[ic], (void*)(*jc), fUsable[ic]?"y":"n");
775   }
776
777   printf("  fYref[0] =%f fYref[1] =%f\n", fYref[0], fYref[1]);
778   printf("  fZref[0] =%f fZref[1] =%f\n", fZref[0], fZref[1]);
779   printf("  fYfit[0] =%f fYfit[1] =%f\n", fYfit[0], fYfit[1]);
780   printf("  fZfit[0] =%f fZfit[1] =%f\n", fZfit[0], fZfit[1]);
781   printf("  fSigmaY =%f\n", fSigmaY);
782   printf("  fSigmaY2=%f\n", fSigmaY2);            
783   printf("  fMeanz  =%f\n", fMeanz);
784   printf("  fZProb  =%f\n", fZProb);
785   printf("  fLabels[0]=%d fLabels[1]=%d\n", fLabels[0], fLabels[1]);
786   printf("  fN      =%d\n", fN);
787   printf("  fN2     =%d (>4 isOK - to be redesigned)\n",fN2);
788   printf("  fNUsed  =%d\n", fNUsed);
789   printf("  fFreq   =%d\n", fFreq);
790   printf("  fNChange=%d\n",  fNChange);
791   printf("  fMPads  =%f\n", fMPads);
792   
793   printf("  fC      =%f\n", fC);        
794   printf("  fCC     =%f\n",fCC);      
795   printf("  fChi2   =%f\n", fChi2);  
796   printf("  fChi2Z  =%f\n", fChi2Z);
797 }
798