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