propagate modifications from time container to the tracklet debug stream
[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(*chamber);
372     ch.SetOwner(); 
373     (*AliTRDtrackerV1::DebugStreamer()) << "AttachClustersIter"
374       << "chamber.="   << &ch
375       << "tracklet.="  << this
376       << "\n";  
377   }
378
379   Float_t  tquality;
380   Double_t kroady = fReconstructor->GetRecoParam() ->GetRoad1y();
381   Double_t kroadz = fPadLength * .5 + 1.;
382   
383   // initialize configuration parameters
384   Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
385   Int_t   niter = kZcorr ? 1 : 2;
386   
387   Double_t yexp, zexp;
388   Int_t ncl = 0;
389   // start seed update
390   for (Int_t iter = 0; iter < niter; iter++) {
391     ncl = 0;
392     for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
393       if(!(layer = chamber->GetTB(iTime))) continue;
394       if(!Int_t(*layer)) continue;
395       
396       // define searching configuration
397       Double_t dxlayer = layer->GetX() - fX0;
398       if(c){
399         zexp = c->GetZ();
400         //Try 2 pad-rows in second iteration
401         if (iter > 0) {
402           zexp = fZref[0] + fZref[1] * dxlayer - zcorr;
403           if (zexp > c->GetZ()) zexp = c->GetZ() + fPadLength*0.5;
404           if (zexp < c->GetZ()) zexp = c->GetZ() - fPadLength*0.5;
405         }
406       } else zexp = fZref[0] + (kZcorr ? fZref[1] * dxlayer : 0.);
407       yexp  = fYref[0] + fYref[1] * dxlayer - zcorr;
408       
409       // Get and register cluster
410       Int_t    index = layer->SearchNearestCluster(yexp, zexp, kroady, kroadz);
411       if (index < 0) continue;
412       AliTRDcluster *cl = (*layer)[index];
413       
414       fIndexes[iTime]  = layer->GetGlobalIndex(index);
415       fClusters[iTime] = cl;
416       fY[iTime]        = cl->GetY();
417       fZ[iTime]        = cl->GetZ();
418       ncl++;
419     }
420     if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d ncl [%d] = %d", iter, fDet, ncl));
421     
422     if(ncl>1){  
423       // calculate length of the time bin (calibration aware)
424       Int_t irp = 0; Float_t x[2]; Int_t tb[2];
425       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
426         if(!fClusters[iTime]) continue;
427         x[irp]  = fClusters[iTime]->GetX();
428         tb[irp] = iTime;
429         irp++;
430         if(irp==2) break;
431       } 
432       fdX = (x[1] - x[0]) / (tb[0] - tb[1]);
433   
434       // update X0 from the clusters (calibration/alignment aware)
435       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
436         if(!(layer = chamber->GetTB(iTime))) continue;
437         if(!layer->IsT0()) continue;
438         if(fClusters[iTime]){ 
439           fX0 = fClusters[iTime]->GetX();
440           break;
441         } else { // we have to infere the position of the anode wire from the other clusters
442           for (Int_t jTime = iTime+1; jTime < AliTRDtrackerV1::GetNTimeBins(); jTime++) {
443             if(!fClusters[jTime]) continue;
444             fX0 = fClusters[jTime]->GetX() + fdX * (jTime - iTime);
445           }
446           break;
447         }
448       } 
449       
450       // update YZ reference point
451       // TODO
452       
453       // update x reference positions (calibration/alignment aware)
454       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
455         if(!fClusters[iTime]) continue;
456         fX[iTime] = fClusters[iTime]->GetX() - fX0;
457       } 
458       
459       AliTRDseed::Update();
460     }
461     if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d nclFit [%d] = %d", iter, fDet, fN2));
462     
463     if(IsOK()){
464       tquality = GetQuality(kZcorr);
465       if(tquality < quality) break;
466       else quality = tquality;
467     }
468     kroadz *= 2.;
469   } // Loop: iter
470   if (!IsOK()) return kFALSE;
471
472   if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=1) CookLabels();
473   UpdateUsed();
474   return kTRUE; 
475 }
476
477 //____________________________________________________________________
478 Bool_t  AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *chamber
479                                       ,Bool_t kZcorr)
480 {
481   //
482   // Projective algorithm to attach clusters to seeding tracklets
483   //
484   // Parameters
485   //
486   // Output
487   //
488   // Detailed description
489   // 1. Collapse x coordinate for the full detector plane
490   // 2. truncated mean on y (r-phi) direction
491   // 3. purge clusters
492   // 4. truncated mean on z direction
493   // 5. purge clusters
494   // 6. fit tracklet
495   //    
496
497   if(!fReconstructor->GetRecoParam() ){
498     AliError("Seed can not be used without a valid RecoParam.");
499     return kFALSE;
500   }
501
502   const Int_t kClusterCandidates = 2 * knTimebins;
503   
504   //define roads
505   Double_t kroady = fReconstructor->GetRecoParam() ->GetRoad1y();
506   Double_t kroadz = fPadLength * 1.5 + 1.;
507   // correction to y for the tilting angle
508   Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
509
510   // working variables
511   AliTRDcluster *clusters[kClusterCandidates];
512   Double_t cond[4], yexp[knTimebins], zexp[knTimebins],
513     yres[kClusterCandidates], zres[kClusterCandidates];
514   Int_t ncl, *index = 0x0, tboundary[knTimebins];
515   
516   // Do cluster projection
517   AliTRDchamberTimeBin *layer = 0x0;
518   Int_t nYclusters = 0; Bool_t kEXIT = kFALSE;
519   for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
520     if(!(layer = chamber->GetTB(iTime))) continue;
521     if(!Int_t(*layer)) continue;
522     
523     fX[iTime] = layer->GetX() - fX0;
524     zexp[iTime] = fZref[0] + fZref[1] * fX[iTime];
525     yexp[iTime] = fYref[0] + fYref[1] * fX[iTime] - zcorr;
526     
527     // build condition and process clusters
528     cond[0] = yexp[iTime] - kroady; cond[1] = yexp[iTime] + kroady;
529     cond[2] = zexp[iTime] - kroadz; cond[3] = zexp[iTime] + kroadz;
530     layer->GetClusters(cond, index, ncl);
531     for(Int_t ic = 0; ic<ncl; ic++){
532       AliTRDcluster *c = layer->GetCluster(index[ic]);
533       clusters[nYclusters] = c;
534       yres[nYclusters++] = c->GetY() - yexp[iTime];
535       if(nYclusters >= kClusterCandidates) {
536         AliWarning(Form("Cluster candidates reached limit %d. Some may be lost.", kClusterCandidates));
537         kEXIT = kTRUE;
538         break;
539       }
540     }
541     tboundary[iTime] = nYclusters;
542     if(kEXIT) break;
543   }
544   
545   // Evaluate truncated mean on the y direction
546   Double_t mean, sigma;
547   AliMathBase::EvaluateUni(nYclusters, yres, mean, sigma, Int_t(nYclusters*.8)-2);
548   // purge cluster candidates
549   Int_t nZclusters = 0;
550   for(Int_t ic = 0; ic<nYclusters; ic++){
551     if(yres[ic] - mean > 4. * sigma){
552       clusters[ic] = 0x0;
553       continue;
554     }
555     zres[nZclusters++] = clusters[ic]->GetZ() - zexp[clusters[ic]->GetLocalTimeBin()];
556   }
557   
558   // Evaluate truncated mean on the z direction
559   AliMathBase::EvaluateUni(nZclusters, zres, mean, sigma, Int_t(nZclusters*.8)-2);
560   // purge cluster candidates
561   for(Int_t ic = 0; ic<nZclusters; ic++){
562     if(zres[ic] - mean > 4. * sigma){
563       clusters[ic] = 0x0;
564       continue;
565     }
566   }
567
568   
569   // Select only one cluster/TimeBin
570   Int_t lastCluster = 0;
571   fN2 = 0;
572   for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
573     ncl = tboundary[iTime] - lastCluster;
574     if(!ncl) continue;
575     Int_t iptr = lastCluster;
576     if(ncl > 1){
577       Float_t dold = 9999.;
578       for(int ic=lastCluster; ic<tboundary[iTime]; ic++){
579         if(!clusters[ic]) continue;
580         Float_t y = yexp[iTime] - clusters[ic]->GetY();
581         Float_t z = zexp[iTime] - clusters[ic]->GetZ();
582         Float_t d = y * y + z * z;
583         if(d > dold) continue;
584         dold = d;
585         iptr = ic;
586       }
587     }
588     fIndexes[iTime]  = chamber->GetTB(iTime)->GetGlobalIndex(iptr);
589     fClusters[iTime] = clusters[iptr];
590     fY[iTime]        = clusters[iptr]->GetY();
591     fZ[iTime]        = clusters[iptr]->GetZ();
592     lastCluster      = tboundary[iTime];
593     fN2++;
594   }
595   
596   // number of minimum numbers of clusters expected for the tracklet
597   Int_t kClmin = Int_t(fReconstructor->GetRecoParam() ->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
598   if (fN2 < kClmin){
599     AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin));
600     fN2 = 0;
601     return kFALSE;
602   }
603
604   // update used clusters
605   fNUsed = 0;
606   for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
607     if(!fClusters[iTime]) continue;
608     if((fClusters[iTime]->IsUsed())) fNUsed++;
609   }
610
611   if (fN2-fNUsed < kClmin){
612     AliWarning(Form("Too many clusters already in use %d (from %d).", fNUsed, fN2));
613     fN2 = 0;
614     return kFALSE;
615   }
616   
617   return kTRUE;
618 }
619
620 //____________________________________________________________________
621 Bool_t AliTRDseedV1::Fit(Bool_t tilt)
622 {
623   //
624   // Linear fit of the tracklet
625   //
626   // Parameters :
627   //
628   // Output :
629   //  True if successful
630   //
631   // Detailed description
632   // 2. Check if tracklet crosses pad row boundary
633   // 1. Calculate residuals in the y (r-phi) direction
634   // 3. Do a Least Square Fit to the data
635   //
636
637   const Int_t kClmin = 8;
638   const Float_t q0 = 100.;
639   const Float_t clSigma0 = 2.E-2;    //[cm]
640   const Float_t clSlopeQ = -1.19E-2; //[1/cm]
641
642   // get track direction
643   Double_t y0   = fYref[0];
644   Double_t dydx = fYref[1]; 
645   Double_t z0   = fZref[0];
646   Double_t dzdx = fZref[1];
647   Double_t yt, zt;
648
649   const Int_t kNtb = AliTRDtrackerV1::GetNTimeBins();
650   AliTRDtrackerV1::AliTRDLeastSquare fitterY, fitterZ;
651
652   // convertion factor from square to gauss distribution for sigma
653   Double_t convert = 1./TMath::Sqrt(12.);
654   
655   // book cluster information
656   Double_t xc[knTimebins], yc[knTimebins], zc[knTimebins], sy[knTimebins], sz[knTimebins];
657   Int_t zRow[knTimebins];
658
659
660   fN = 0;
661   AliTRDcluster *c=0x0, **jc = &fClusters[0];
662   for (Int_t ic=0; ic<kNtb; ic++, ++jc) {
663     zRow[ic] = -1;
664     xc[ic]  = -1.;
665     yc[ic]  = 999.;
666     zc[ic]  = 999.;
667     sy[ic]  = 0.;
668     sz[ic]  = 0.;
669     if(!(c = (*jc))) continue;
670     if(!c->IsInChamber()) continue;
671     Float_t w = 1.;
672     if(c->GetNPads()>4) w = .5;
673     if(c->GetNPads()>5) w = .2;
674     zRow[fN] = c->GetPadRow();
675     xc[fN]   = fX0 - c->GetX();
676     yc[fN]   = c->GetY();
677     zc[fN]   = c->GetZ();
678
679     // extrapolated y value for the track
680     yt = y0 - xc[fN]*dydx; 
681     // extrapolated z value for the track
682     zt = z0 - xc[fN]*dzdx; 
683     // tilt correction
684     if(tilt) yc[fN] -= fTilt*(zc[fN] - zt); 
685
686     // elaborate cluster error
687     Float_t qr = c->GetQ() - q0;
688     sy[fN]   = qr < 0. ? clSigma0*TMath::Exp(clSlopeQ*qr) : clSigma0;
689
690     fitterY.AddPoint(&xc[fN], yc[fN]-yt, sy[fN]);
691
692     sz[fN]   = fPadLength*convert;
693     fitterZ.AddPoint(&xc[fN], zc[fN], sz[fN]);
694     fN++;
695   }
696   // to few clusters
697   if (fN < kClmin) return kFALSE; 
698
699   // fit XY plane
700   fitterY.Eval();
701   fYfit[0] = y0+fitterY.GetFunctionParameter(0);
702   fYfit[1] = dydx-fitterY.GetFunctionParameter(1);
703
704   // check par row crossing
705   Int_t zN[2*AliTRDseed::knTimebins];
706   Int_t nz = AliTRDtrackerV1::Freq(fN, zRow, zN, kFALSE);
707   // more than one pad row crossing
708   if(nz>2) return kFALSE; 
709
710
711   // determine z offset of the fit
712   Float_t zslope = 0.;
713   Int_t nchanges = 0, nCross = 0;
714   if(nz==2){ // tracklet is crossing pad row
715     // Find the break time allowing one chage on pad-rows
716     // with maximal number of accepted clusters
717     Int_t padRef = zRow[0];
718     for (Int_t ic=1; ic<fN; ic++) {
719       if(zRow[ic] == padRef) continue;
720       
721       // debug
722       if(zRow[ic-1] == zRow[ic]){
723         printf("ERROR in pad row change!!!\n");
724       }
725     
726       // evaluate parameters of the crossing point
727       Float_t sx = (xc[ic-1] - xc[ic])*convert;
728       fCross[0] = .5 * (xc[ic-1] + xc[ic]);
729       fCross[2] = .5 * (zc[ic-1] + zc[ic]);
730       fCross[3] = TMath::Max(dzdx * sx, .01);
731       zslope    = zc[ic-1] > zc[ic] ? 1. : -1.;
732       padRef    = zRow[ic];
733       nCross    = ic;
734       nchanges++;
735     }
736   }
737
738   // condition on nCross and reset nchanges TODO
739
740   if(nchanges==1){
741     if(dzdx * zslope < 0.){
742       AliInfo("tracklet direction does not correspond to the track direction. TODO.");
743     }
744     SetBit(kRowCross, kTRUE); // mark pad row crossing
745     fitterZ.AddPoint(&fCross[0], fCross[2], fCross[3]);
746     fitterZ.Eval();
747     //zc[nc] = fitterZ.GetFunctionParameter(0); 
748     fCross[1] = fYfit[0] - fCross[0] * fYfit[1];
749     fCross[0] = fX0 - fCross[0];
750   } else if(nchanges > 1){ // debug
751     AliError("N pad row crossing > 1.");
752     return kFALSE;
753   }
754
755   UpdateUsed();
756
757   return kTRUE;
758 }
759
760
761 //___________________________________________________________________
762 void AliTRDseedV1::Print(Option_t*) const
763 {
764   //
765   // Printing the seedstatus
766   //
767
768   AliInfo(Form("Tracklet X0[%7.2f] Det[%d]", fX0, fDet));
769   printf("  Tilt[%+6.2f] PadLength[%5.2f]\n", fTilt, fPadLength);
770   AliTRDcluster* const* jc = &fClusters[0];
771   for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++, jc++) {
772     if(!(*jc)) continue;
773     printf("  %2d X[%7.2f] Y[%7.2f] Z[%7.2f] Idx[%d] c[%p] usable[%s]\n", 
774       ic, (*jc)->GetX(), (*jc)->GetY(), (*jc)->GetZ(), 
775       fIndexes[ic], (void*)(*jc), fUsable[ic]?"y":"n");
776   }
777
778   printf("  fYref[0] =%f fYref[1] =%f\n", fYref[0], fYref[1]);
779   printf("  fZref[0] =%f fZref[1] =%f\n", fZref[0], fZref[1]);
780   printf("  fYfit[0] =%f fYfit[1] =%f\n", fYfit[0], fYfit[1]);
781   printf("  fZfit[0] =%f fZfit[1] =%f\n", fZfit[0], fZfit[1]);
782   printf("  fSigmaY =%f\n", fSigmaY);
783   printf("  fSigmaY2=%f\n", fSigmaY2);            
784   printf("  fMeanz  =%f\n", fMeanz);
785   printf("  fZProb  =%f\n", fZProb);
786   printf("  fLabels[0]=%d fLabels[1]=%d\n", fLabels[0], fLabels[1]);
787   printf("  fN      =%d\n", fN);
788   printf("  fN2     =%d (>4 isOK - to be redesigned)\n",fN2);
789   printf("  fNUsed  =%d\n", fNUsed);
790   printf("  fFreq   =%d\n", fFreq);
791   printf("  fNChange=%d\n",  fNChange);
792   printf("  fMPads  =%f\n", fMPads);
793   
794   printf("  fC      =%f\n", fC);        
795   printf("  fCC     =%f\n",fCC);      
796   printf("  fChi2   =%f\n", fChi2);  
797   printf("  fChi2Z  =%f\n", fChi2Z);
798 }
799