caa076758c81bbe04fd693e9a65d19a1b5c29626
[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 #include "AliCDBManager.h"
36 #include "AliTracker.h"
37
38 #include "AliTRDpadPlane.h"
39 #include "AliTRDcluster.h"
40 #include "AliTRDseedV1.h"
41 #include "AliTRDtrackV1.h"
42 #include "AliTRDcalibDB.h"
43 #include "AliTRDchamberTimeBin.h"
44 #include "AliTRDtrackingChamber.h"
45 #include "AliTRDtrackerV1.h"
46 #include "AliTRDReconstructor.h"
47 #include "AliTRDrecoParam.h"
48 #include "AliTRDCommonParam.h"
49
50 #include "Cal/AliTRDCalPID.h"
51 #include "Cal/AliTRDCalROC.h"
52 #include "Cal/AliTRDCalDet.h"
53
54 ClassImp(AliTRDseedV1)
55
56 //____________________________________________________________________
57 AliTRDseedV1::AliTRDseedV1(Int_t det) 
58   :AliTRDseed()
59   ,fReconstructor(0x0)
60   ,fClusterIter(0x0)
61   ,fClusterIdx(0)
62   ,fDet(det)
63   ,fMom(0.)
64   ,fSnp(0.)
65   ,fTgl(0.)
66   ,fdX(0.)
67   ,fXref(0.)
68   ,fExB(0.)
69 {
70   //
71   // Constructor
72   //
73   for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = 0.;
74   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec]  = -1.;
75   fRefCov[0] = 1.; fRefCov[1] = 0.; fRefCov[2] = 1.;
76   // covariance matrix [diagonal]
77   // default sy = 200um and sz = 2.3 cm 
78   fCov[0] = 4.e-4; fCov[1] = 0.; fCov[2] = 5.3; 
79 }
80
81 //____________________________________________________________________
82 AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
83   :AliTRDseed((AliTRDseed&)ref)
84   ,fReconstructor(ref.fReconstructor)
85   ,fClusterIter(0x0)
86   ,fClusterIdx(0)
87   ,fDet(ref.fDet)
88   ,fMom(ref.fMom)
89   ,fSnp(ref.fSnp)
90   ,fTgl(ref.fTgl)
91   ,fdX(ref.fdX)
92   ,fXref(ref.fXref)
93   ,fExB(ref.fExB)
94 {
95   //
96   // Copy Constructor performing a deep copy
97   //
98
99   //printf("AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &)\n");
100   SetBit(kOwner, kFALSE);
101   for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = ref.fdEdx[islice];
102   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = ref.fProb[ispec];
103   memcpy(fRefCov, ref.fRefCov, 3*sizeof(Double_t));
104   memcpy(fCov, ref.fCov, 3*sizeof(Double_t));
105 }
106
107
108 //____________________________________________________________________
109 AliTRDseedV1& AliTRDseedV1::operator=(const AliTRDseedV1 &ref)
110 {
111   //
112   // Assignment Operator using the copy function
113   //
114
115   if(this != &ref){
116     ref.Copy(*this);
117   }
118   SetBit(kOwner, kFALSE);
119
120   return *this;
121
122 }
123
124 //____________________________________________________________________
125 AliTRDseedV1::~AliTRDseedV1()
126 {
127   //
128   // Destructor. The RecoParam object belongs to the underlying tracker.
129   //
130
131   //printf("I-AliTRDseedV1::~AliTRDseedV1() : Owner[%s]\n", IsOwner()?"YES":"NO");
132
133   if(IsOwner()) 
134     for(int itb=0; itb<knTimebins; itb++){
135       if(!fClusters[itb]) continue; 
136       //AliInfo(Form("deleting c %p @ %d", fClusters[itb], itb));
137       delete fClusters[itb];
138       fClusters[itb] = 0x0;
139     }
140 }
141
142 //____________________________________________________________________
143 void AliTRDseedV1::Copy(TObject &ref) const
144 {
145   //
146   // Copy function
147   //
148
149   //AliInfo("");
150   AliTRDseedV1 &target = (AliTRDseedV1 &)ref; 
151
152   target.fClusterIter   = 0x0;
153   target.fClusterIdx    = 0;
154   target.fDet           = fDet;
155   target.fMom           = fMom;
156   target.fSnp           = fSnp;
157   target.fTgl           = fTgl;
158   target.fdX            = fdX;
159   target.fXref          = fXref;
160   target.fExB           = fExB;
161   target.fReconstructor = fReconstructor;
162   
163   for(int islice=0; islice < knSlices; islice++) target.fdEdx[islice] = fdEdx[islice];
164   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) target.fProb[ispec] = fProb[ispec];
165   memcpy(target.fRefCov, fRefCov, 3*sizeof(Double_t));
166   memcpy(target.fCov, fCov, 3*sizeof(Double_t));
167   
168   AliTRDseed::Copy(target);
169 }
170
171
172 //____________________________________________________________
173 Bool_t AliTRDseedV1::Init(AliTRDtrackV1 *track)
174 {
175 // Initialize this tracklet using the track information
176 //
177 // Parameters:
178 //   track - the TRD track used to initialize the tracklet
179 // 
180 // Detailed description
181 // The function sets the starting point and direction of the
182 // tracklet according to the information from the TRD track.
183 // 
184 // Caution
185 // The TRD track has to be propagated to the beginning of the
186 // chamber where the tracklet will be constructed
187 //
188
189   Double_t y, z; 
190   if(!track->GetProlongation(fX0, y, z)) return kFALSE;
191   UpDate(track);
192   return kTRUE;
193 }
194
195
196 //____________________________________________________________________
197 void AliTRDseedV1::UpDate(const AliTRDtrackV1 *trk)
198
199   // update tracklet reference position from the TRD track
200   // Funny name to avoid the clash with the function AliTRDseed::Update() (has to be made obselete)
201
202   fSnp = trk->GetSnp();
203   fTgl = trk->GetTgl();
204   fMom = trk->GetP();
205   fYref[1] = fSnp/(1. - fSnp*fSnp);
206   fZref[1] = fTgl;
207   SetCovRef(trk->GetCovariance());
208
209   Double_t dx = trk->GetX() - fX0;
210   fYref[0] = trk->GetY() - dx*fYref[1];
211   fZref[0] = trk->GetZ() - dx*fZref[1];
212 }
213
214 //____________________________________________________________________
215 void AliTRDseedV1::CookdEdx(Int_t nslices)
216 {
217 // Calculates average dE/dx for all slices and store them in the internal array fdEdx. 
218 //
219 // Parameters:
220 //  nslices : number of slices for which dE/dx should be calculated
221 // Output:
222 //  store results in the internal array fdEdx. This can be accessed with the method
223 //  AliTRDseedV1::GetdEdx()
224 //
225 // Detailed description
226 // Calculates average dE/dx for all slices. Depending on the PID methode 
227 // the number of slices can be 3 (LQ) or 8(NN). 
228 // The calculation of dQ/dl are done using the tracklet fit results (see AliTRDseedV1::GetdQdl(Int_t))
229 //
230 // The following effects are included in the calculation:
231 // 1. calibration values for t0 and vdrift (using x coordinate to calculate slice)
232 // 2. cluster sharing (optional see AliTRDrecoParam::SetClusterSharing())
233 // 3. cluster size
234 //
235
236   Int_t nclusters[knSlices];
237   for(int i=0; i<knSlices; i++){ 
238     fdEdx[i]     = 0.;
239     nclusters[i] = 0;
240   }
241   const Double_t kDriftLength = (.5 * AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
242
243   AliTRDcluster *c = 0x0;
244   for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++){
245     if(!(c = fClusters[ic]) && !(c = fClusters[ic+kNtb])) continue;
246     Float_t dx = TMath::Abs(fX0 - c->GetX());
247     
248     // Filter clusters for dE/dx calculation
249     
250     // 1.consider calibration effects for slice determination
251     Int_t slice;
252     if(dx<kDriftLength){ // TODO should be replaced by c->IsInChamber() 
253       slice = Int_t(dx * nslices / kDriftLength);
254     } else slice = c->GetX() < fX0 ? nslices-1 : 0;
255
256
257     // 2. take sharing into account
258     Float_t w = c->IsShared() ? .5 : 1.;
259     
260     // 3. take into account large clusters TODO
261     //w *= c->GetNPads() > 3 ? .8 : 1.;
262     
263     //CHECK !!!
264     fdEdx[slice]   += w * GetdQdl(ic); //fdQdl[ic];
265     nclusters[slice]++;
266   } // End of loop over clusters
267
268   //if(fReconstructor->GetPIDMethod() == AliTRDReconstructor::kLQPID){
269   if(nslices == AliTRDpidUtil::kLQslices){
270   // calculate mean charge per slice (only LQ PID)
271     for(int is=0; is<nslices; is++){ 
272       if(nclusters[is]) fdEdx[is] /= nclusters[is];
273     }
274   }
275 }
276
277 //____________________________________________________________________
278 void AliTRDseedV1::GetClusterXY(const AliTRDcluster *c, Double_t &x, Double_t &y)
279 {
280 // Return corrected position of the cluster taking into 
281 // account variation of the drift velocity with drift length.
282
283
284   // drift velocity correction TODO to be moved to the clusterizer
285   const Float_t cx[] = {
286     -9.6280e-02, 1.3091e-01,-1.7415e-02,-9.9221e-02,-1.2040e-01,-9.5493e-02,
287     -5.0041e-02,-1.6726e-02, 3.5756e-03, 1.8611e-02, 2.6378e-02, 3.3823e-02,
288      3.4811e-02, 3.5282e-02, 3.5386e-02, 3.6047e-02, 3.5201e-02, 3.4384e-02,
289      3.2864e-02, 3.1932e-02, 3.2051e-02, 2.2539e-02,-2.5154e-02,-1.2050e-01,
290     -1.2050e-01
291   };
292
293   // PRF correction TODO to be replaced by the gaussian 
294   // approximation with full error parametrization and // moved to the clusterizer
295   const Float_t cy[AliTRDgeometry::kNlayer][3] = {
296     { 4.014e-04, 8.605e-03, -6.880e+00},
297     {-3.061e-04, 9.663e-03, -6.789e+00},
298     { 1.124e-03, 1.105e-02, -6.825e+00},
299     {-1.527e-03, 1.231e-02, -6.777e+00},
300     { 2.150e-03, 1.387e-02, -6.783e+00},
301     {-1.296e-03, 1.486e-02, -6.825e+00}
302   }; 
303
304   Int_t ily = AliTRDgeometry::GetLayer(c->GetDetector());
305   x = c->GetX() - cx[c->GetLocalTimeBin()];
306   y = c->GetY() + cy[ily][0] + cy[ily][1] * TMath::Sin(cy[ily][2] * c->GetCenter());
307   return;
308 }
309
310 //____________________________________________________________________
311 Float_t AliTRDseedV1::GetdQdl(Int_t ic) const
312 {
313 // Using the linear approximation of the track inside one TRD chamber (TRD tracklet) 
314 // the charge per unit length can be written as:
315 // BEGIN_LATEX
316 // #frac{dq}{dl} = #frac{q_{c}}{dx * #sqrt{1 + #(){#frac{dy}{dx}}^{2}_{fit} + #(){#frac{dy}{dx}}^{2}_{ref}}}
317 // END_LATEX
318 // where qc is the total charge collected in the current time bin and dx is the length 
319 // of the time bin. For the moment (Jan 20 2009) only pad row cross corrections are 
320 // considered for the charge but none are applied for drift velocity variations along 
321 // the drift region or assymetry of the TRF
322 // 
323 // Author : Alex Bercuci <A.Bercuci@gsi.de>
324 //
325   Float_t dq = 0.;
326   if(fClusters[ic]) dq += TMath::Abs(fClusters[ic]->GetQ());
327   if(fClusters[ic+kNtb]) dq += TMath::Abs(fClusters[ic+kNtb]->GetQ());
328   if(dq<1.e-3 || fdX < 1.e-3) return 0.;
329
330   return dq/fdX/TMath::Sqrt(1. + fYfit[1]*fYfit[1] + fZref[1]*fZref[1]);
331 }
332
333 //____________________________________________________________________
334 Double_t* AliTRDseedV1::GetProbability()
335 {       
336 // Fill probability array for tracklet from the DB.
337 //
338 // Parameters
339 //
340 // Output
341 //   returns pointer to the probability array and 0x0 if missing DB access 
342 //
343 // Detailed description
344
345   
346   // retrive calibration db
347   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
348   if (!calibration) {
349     AliError("No access to calibration data");
350     return 0x0;
351   }
352
353   if (!fReconstructor) {
354     AliError("Reconstructor not set.");
355     return 0x0;
356   }
357
358   // Retrieve the CDB container class with the parametric detector response
359   const AliTRDCalPID *pd = calibration->GetPIDObject(fReconstructor->GetPIDMethod());
360   if (!pd) {
361     AliError("No access to AliTRDCalPID object");
362     return 0x0;
363   }
364   //AliInfo(Form("Method[%d] : %s", fReconstructor->GetRecoParam() ->GetPIDMethod(), pd->IsA()->GetName()));
365
366   // calculate tracklet length TO DO
367   Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
368   /// TMath::Sqrt((1.0 - fSnp[iPlane]*fSnp[iPlane]) / (1.0 + fTgl[iPlane]*fTgl[iPlane]));
369   
370   //calculate dE/dx
371   CookdEdx(fReconstructor->GetNdEdxSlices());
372   
373   // Sets the a priori probabilities
374   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
375     fProb[ispec] = pd->GetProbability(ispec, fMom, &fdEdx[0], length, GetPlane());      
376   }
377
378   return &fProb[0];
379 }
380
381 //____________________________________________________________________
382 Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
383 {
384   //
385   // Returns a quality measurement of the current seed
386   //
387
388   Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
389   return 
390       .5 * TMath::Abs(18.0 - fN2)
391     + 10.* TMath::Abs(fYfit[1] - fYref[1])
392     + 5. * TMath::Abs(fYfit[0] - fYref[0] + zcorr)
393     + 2. * TMath::Abs(fMeanz - fZref[0]) / fPadLength;
394 }
395
396 //____________________________________________________________________
397 void AliTRDseedV1::GetCovAt(Double_t x, Double_t *cov) const
398 {
399 // Computes covariance in the y-z plane at radial point x (in tracking coordinates) 
400 // and returns the results in the preallocated array cov[3] as :
401 //   cov[0] = Var(y)
402 //   cov[1] = Cov(yz)
403 //   cov[2] = Var(z)
404 //
405 // Details
406 //
407 // For the linear transformation
408 // BEGIN_LATEX
409 // Y = T_{x} X^{T}
410 // END_LATEX
411 //   The error propagation has the general form
412 // BEGIN_LATEX
413 // C_{Y} = T_{x} C_{X} T_{x}^{T} 
414 // END_LATEX
415 //  We apply this formula 2 times. First to calculate the covariance of the tracklet 
416 // at point x we consider: 
417 // BEGIN_LATEX
418 // T_{x} = (1 x); X=(y0 dy/dx); C_{X}=#(){#splitline{Var(y0) Cov(y0, dy/dx)}{Cov(y0, dy/dx) Var(dy/dx)}} 
419 // END_LATEX
420 // and secondly to take into account the tilt angle
421 // BEGIN_LATEX
422 // T_{#alpha} = #(){#splitline{cos(#alpha) __ sin(#alpha)}{-sin(#alpha) __ cos(#alpha)}}; X=(y z); C_{X}=#(){#splitline{Var(y)    0}{0   Var(z)}} 
423 // END_LATEX
424 //
425 // using simple trigonometrics one can write for this last case
426 // BEGIN_LATEX
427 // C_{Y}=#frac{1}{1+tg^{2}#alpha} #(){#splitline{(#sigma_{y}^{2}+tg^{2}#alpha#sigma_{z}^{2}) __ tg#alpha(#sigma_{z}^{2}-#sigma_{y}^{2})}{tg#alpha(#sigma_{z}^{2}-#sigma_{y}^{2}) __ (#sigma_{z}^{2}+tg^{2}#alpha#sigma_{y}^{2})}} 
428 // END_LATEX
429 // which can be aproximated for small alphas (2 deg) with
430 // BEGIN_LATEX
431 // C_{Y}=#(){#splitline{#sigma_{y}^{2} __ (#sigma_{z}^{2}-#sigma_{y}^{2})tg#alpha}{((#sigma_{z}^{2}-#sigma_{y}^{2})tg#alpha __ #sigma_{z}^{2}}} 
432 // END_LATEX
433 //
434 // before applying the tilt rotation we also apply systematic uncertainties to the tracklet 
435 // position which can be tunned from outside via the AliTRDrecoParam::SetSysCovMatrix(). They might 
436 // account for extra misalignment/miscalibration uncertainties. 
437 //
438 // Author :
439 // Alex Bercuci <A.Bercuci@gsi.de> 
440 // Date : Jan 8th 2009
441 //
442
443
444   Double_t xr     = fX0-x; 
445   Double_t sy2    = fCov[0] +2.*xr*fCov[1] + xr*xr*fCov[2];
446   Double_t sz2    = fPadLength*fPadLength/12.;
447
448   // insert systematic uncertainties
449   Double_t sys[15];
450   fReconstructor->GetRecoParam()->GetSysCovMatrix(sys);
451   sy2 += sys[0];
452   sz2 += sys[1];
453
454   // rotate covariance matrix
455   Double_t t2 = fTilt*fTilt;
456   Double_t correction = 1./(1. + t2);
457   cov[0] = (sy2+t2*sz2)*correction;
458   cov[1] = fTilt*(sz2 - sy2)*correction;
459   cov[2] = (t2*sy2+sz2)*correction;
460 }
461
462
463 //____________________________________________________________________
464 void AliTRDseedV1::SetExB()
465 {
466 // Retrive the tg(a_L) from OCDB. The following information are used
467 //  - detector index
468 //  - column and row position of first attached cluster. 
469 // 
470 // If no clusters are attached to the tracklet a random central chamber 
471 // position (c=70, r=7) will be used to retrieve the drift velocity.
472 //
473 // Author :
474 // Alex Bercuci <A.Bercuci@gsi.de> 
475 // Date : Jan 8th 2009
476 //
477
478   AliCDBManager *cdb = AliCDBManager::Instance();
479   if(cdb->GetRun() < 0){
480     AliError("OCDB manager not properly initialized");
481     return;
482   }
483
484   AliTRDcalibDB *fCalib = AliTRDcalibDB::Instance();
485   AliTRDCalROC  *fCalVdriftROC = fCalib->GetVdriftROC(fDet);
486   const AliTRDCalDet  *fCalVdriftDet = fCalib->GetVdriftDet();
487
488   Int_t col = 70, row = 7;
489   AliTRDcluster **c = &fClusters[0];
490   if(fN){ 
491     Int_t ic = 0;
492     while (ic<AliTRDseed::knTimebins && !(*c)){ic++; c++;} 
493     if(*c){
494       col = (*c)->GetPadCol();
495       row = (*c)->GetPadRow();
496     }
497   }
498
499   Double_t vd = fCalVdriftDet->GetValue(fDet) * fCalVdriftROC->GetValue(col, row);
500   fExB   = AliTRDCommonParam::Instance()->GetOmegaTau(vd);
501 }
502
503 //____________________________________________________________________
504 void AliTRDseedV1::SetOwner()
505 {
506   //AliInfo(Form("own [%s] fOwner[%s]", own?"YES":"NO", fOwner?"YES":"NO"));
507   
508   if(TestBit(kOwner)) return;
509   for(int ic=0; ic<knTimebins; ic++){
510     if(!fClusters[ic]) continue;
511     fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
512   }
513   SetBit(kOwner);
514 }
515
516 //____________________________________________________________________
517 Bool_t  AliTRDseedV1::AttachClustersIter(AliTRDtrackingChamber *chamber, Float_t quality, Bool_t kZcorr, AliTRDcluster *c)
518 {
519   //
520   // Iterative process to register clusters to the seed.
521   // In iteration 0 we try only one pad-row and if quality not
522   // sufficient we try 2 pad-rows (about 5% of tracks cross 2 pad-rows)
523   //
524   // debug level 7
525   //
526   
527   if(!fReconstructor->GetRecoParam() ){
528     AliError("Seed can not be used without a valid RecoParam.");
529     return kFALSE;
530   }
531
532   AliTRDchamberTimeBin *layer = 0x0;
533   if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7){
534     AliTRDtrackingChamber ch(*chamber);
535     ch.SetOwner(); 
536     TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker);
537     cstreamer << "AttachClustersIter"
538       << "chamber.="   << &ch
539       << "tracklet.="  << this
540       << "\n";  
541   }
542
543   Float_t  tquality;
544   Double_t kroady = fReconstructor->GetRecoParam() ->GetRoad1y();
545   Double_t kroadz = fPadLength * .5 + 1.;
546   
547   // initialize configuration parameters
548   Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
549   Int_t   niter = kZcorr ? 1 : 2;
550   
551   Double_t yexp, zexp;
552   Int_t ncl = 0;
553   // start seed update
554   for (Int_t iter = 0; iter < niter; iter++) {
555     ncl = 0;
556     for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
557       if(!(layer = chamber->GetTB(iTime))) continue;
558       if(!Int_t(*layer)) continue;
559       
560       // define searching configuration
561       Double_t dxlayer = layer->GetX() - fX0;
562       if(c){
563         zexp = c->GetZ();
564         //Try 2 pad-rows in second iteration
565         if (iter > 0) {
566           zexp = fZref[0] + fZref[1] * dxlayer - zcorr;
567           if (zexp > c->GetZ()) zexp = c->GetZ() + fPadLength*0.5;
568           if (zexp < c->GetZ()) zexp = c->GetZ() - fPadLength*0.5;
569         }
570       } else zexp = fZref[0] + (kZcorr ? fZref[1] * dxlayer : 0.);
571       yexp  = fYref[0] + fYref[1] * dxlayer - zcorr;
572       
573       // Get and register cluster
574       Int_t    index = layer->SearchNearestCluster(yexp, zexp, kroady, kroadz);
575       if (index < 0) continue;
576       AliTRDcluster *cl = (*layer)[index];
577       
578       fIndexes[iTime]  = layer->GetGlobalIndex(index);
579       fClusters[iTime] = cl;
580       fY[iTime]        = cl->GetY();
581       fZ[iTime]        = cl->GetZ();
582       ncl++;
583     }
584     if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d ncl [%d] = %d", iter, fDet, ncl));
585     
586     if(ncl>1){  
587       // calculate length of the time bin (calibration aware)
588       Int_t irp = 0; Float_t x[2]; Int_t tb[2];
589       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
590         if(!fClusters[iTime]) continue;
591         x[irp]  = fClusters[iTime]->GetX();
592         tb[irp] = iTime;
593         irp++;
594         if(irp==2) break;
595       } 
596       fdX = (x[1] - x[0]) / (tb[0] - tb[1]);
597   
598       // update X0 from the clusters (calibration/alignment aware)
599       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
600         if(!(layer = chamber->GetTB(iTime))) continue;
601         if(!layer->IsT0()) continue;
602         if(fClusters[iTime]){ 
603           fX0 = fClusters[iTime]->GetX();
604           break;
605         } else { // we have to infere the position of the anode wire from the other clusters
606           for (Int_t jTime = iTime+1; jTime < AliTRDtrackerV1::GetNTimeBins(); jTime++) {
607             if(!fClusters[jTime]) continue;
608             fX0 = fClusters[jTime]->GetX() + fdX * (jTime - iTime);
609             break;
610           }
611         }
612       } 
613       
614       // update YZ reference point
615       // TODO
616       
617       // update x reference positions (calibration/alignment aware)
618       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
619         if(!fClusters[iTime]) continue;
620         fX[iTime] = fX0 - fClusters[iTime]->GetX();
621       } 
622       
623       AliTRDseed::Update();
624     }
625     if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d nclFit [%d] = %d", iter, fDet, fN2));
626     
627     if(IsOK()){
628       tquality = GetQuality(kZcorr);
629       if(tquality < quality) break;
630       else quality = tquality;
631     }
632     kroadz *= 2.;
633   } // Loop: iter
634   if (!IsOK()) return kFALSE;
635
636   if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=1) CookLabels();
637
638   // set ExB angle
639   SetExB();
640   UpdateUsed();
641   return kTRUE; 
642 }
643
644 //____________________________________________________________________
645 Bool_t  AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *chamber, Bool_t tilt)
646 {
647   //
648   // Projective algorithm to attach clusters to seeding tracklets
649   //
650   // Parameters
651   //
652   // Output
653   //
654   // Detailed description
655   // 1. Collapse x coordinate for the full detector plane
656   // 2. truncated mean on y (r-phi) direction
657   // 3. purge clusters
658   // 4. truncated mean on z direction
659   // 5. purge clusters
660   // 6. fit tracklet
661   //    
662   Bool_t kPRINT = kFALSE;
663   if(!fReconstructor->GetRecoParam() ){
664     AliError("Seed can not be used without a valid RecoParam.");
665     return kFALSE;
666   }
667   // Initialize reco params for this tracklet
668   // 1. first time bin in the drift region
669   Int_t t0 = 4;
670   Int_t kClmin = Int_t(fReconstructor->GetRecoParam() ->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
671
672   Double_t syRef  = TMath::Sqrt(fRefCov[0]);
673   //define roads
674   Double_t kroady = 1.; 
675   //fReconstructor->GetRecoParam() ->GetRoad1y();
676   Double_t kroadz = fPadLength * 1.5 + 1.;
677   if(kPRINT) printf("AttachClusters() sy[%f] road[%f]\n", syRef, kroady);
678
679   // working variables
680   const Int_t kNrows = 16;
681   AliTRDcluster *clst[kNrows][knTimebins];
682   Double_t cond[4], dx, dy, yt, zt,
683     yres[kNrows][knTimebins];
684   Int_t idxs[kNrows][knTimebins], ncl[kNrows], ncls = 0;
685   memset(ncl, 0, kNrows*sizeof(Int_t));
686   memset(clst, 0, kNrows*knTimebins*sizeof(AliTRDcluster*));
687
688   // Do cluster projection
689   AliTRDcluster *c = 0x0;
690   AliTRDchamberTimeBin *layer = 0x0;
691   Bool_t kBUFFER = kFALSE;
692   for (Int_t it = 0; it < AliTRDtrackerV1::GetNTimeBins(); it++) {
693     if(!(layer = chamber->GetTB(it))) continue;
694     if(!Int_t(*layer)) continue;
695     
696     dx   = fX0 - layer->GetX();
697     yt = fYref[0] - fYref[1] * dx;
698     zt = fZref[0] - fZref[1] * dx;
699     if(kPRINT) printf("\t%2d dx[%f] yt[%f] zt[%f]\n", it, dx, yt, zt);
700
701     // select clusters on a 5 sigmaKalman level
702     cond[0] = yt; cond[2] = kroady;
703     cond[1] = zt; cond[3] = kroadz;
704     Int_t n=0, idx[6];
705     layer->GetClusters(cond, idx, n, 6);
706     for(Int_t ic = n; ic--;){
707       c  = (*layer)[idx[ic]];
708       dy = yt - c->GetY();
709       dy += tilt ? fTilt * (c->GetZ() - zt) : 0.;
710       // select clusters on a 3 sigmaKalman level
711 /*      if(tilt && TMath::Abs(dy) > 3.*syRef){ 
712         printf("too large !!!\n");
713         continue;
714       }*/
715       Int_t r = c->GetPadRow();
716       if(kPRINT) printf("\t\t%d dy[%f] yc[%f] r[%d]\n", ic, TMath::Abs(dy), c->GetY(), r);
717       clst[r][ncl[r]] = c;
718       idxs[r][ncl[r]] = idx[ic];
719       yres[r][ncl[r]] = dy;
720       ncl[r]++; ncls++;
721
722       if(ncl[r] >= knTimebins) {
723         AliWarning(Form("Cluster candidates reached limit %d. Some may be lost.", knTimebins));
724         kBUFFER = kTRUE;
725         break;
726       }
727     }
728     if(kBUFFER) break;
729   }
730   if(kPRINT) printf("Found %d clusters\n", ncls);
731   if(ncls<kClmin) return kFALSE;
732  
733   // analyze each row individualy
734   Double_t mean, syDis;
735   Int_t nrow[] = {0, 0, 0}, nr = 0, lr=-1;
736   for(Int_t ir=kNrows; ir--;){
737     if(!(ncl[ir])) continue;
738     if(lr>0 && lr-ir != 1){
739       if(kPRINT) printf("W - gap in rows attached !!\n"); 
740     }
741     if(kPRINT) printf("\tir[%d] lr[%d] n[%d]\n", ir, lr, ncl[ir]);
742     // Evaluate truncated mean on the y direction
743     if(ncl[ir] > 3) AliMathBase::EvaluateUni(ncl[ir], yres[ir], mean, syDis, Int_t(ncl[ir]*.8));
744     else {
745       mean = 0.; syDis = 0.;
746     } 
747
748     // TODO check mean and sigma agains cluster resolution !!
749     if(kPRINT) printf("\tr[%2d] m[%f %5.3fsigma] s[%f]\n", ir, mean, TMath::Abs(mean/syRef), syDis);
750     // select clusters on a 3 sigmaDistr level
751     Bool_t kFOUND = kFALSE;
752     for(Int_t ic = ncl[ir]; ic--;){
753       if(yres[ir][ic] - mean > 3. * syDis){ 
754         clst[ir][ic] = 0x0; continue;
755       }
756       nrow[nr]++; kFOUND = kTRUE;
757     }
758     // exit loop
759     if(kFOUND) nr++; 
760     lr = ir; if(nr>=3) break;
761   }
762   if(kPRINT) printf("lr[%d] nr[%d] nrow[0]=%d nrow[1]=%d nrow[2]=%d\n", lr, nr, nrow[0], nrow[1], nrow[2]);
763
764   // classify cluster rows
765   Int_t row = -1;
766   switch(nr){
767   case 1:
768     row = lr;
769     break;
770   case 2:
771     SetBit(kRowCross, kTRUE); // mark pad row crossing
772     if(nrow[0] > nrow[1]){ row = lr+1; lr = -1;}
773     else{ 
774       row = lr; lr = 1;
775       nrow[2] = nrow[1];
776       nrow[1] = nrow[0];
777       nrow[0] = nrow[2];
778     }
779     break;
780   case 3:
781     SetBit(kRowCross, kTRUE); // mark pad row crossing
782     break;
783   }
784   if(kPRINT) printf("\trow[%d] n[%d]\n\n", row, nrow[0]);
785   if(row<0) return kFALSE;
786
787   // Select and store clusters 
788   // We should consider here :
789   //  1. How far is the chamber boundary
790   //  2. How big is the mean
791   fN2 = 0;
792   for (Int_t ir = 0; ir < nr; ir++) {
793     Int_t jr = row + ir*lr; 
794     if(kPRINT) printf("\tattach %d clusters for row %d\n", ncl[jr], jr);
795     for (Int_t ic = 0; ic < ncl[jr]; ic++) {
796       if(!(c = clst[jr][ic])) continue;
797       Int_t it = c->GetPadTime();
798       // TODO proper indexing of clusters !!
799       fIndexes[it+35*ir]  = chamber->GetTB(it)->GetGlobalIndex(idxs[jr][ic]);
800       fClusters[it+35*ir] = c;
801   
802       //printf("\tid[%2d] it[%d] idx[%d]\n", ic, it, fIndexes[it]);
803   
804       fN2++;
805     }
806   }  
807
808   // number of minimum numbers of clusters expected for the tracklet
809   if (fN2 < kClmin){
810     AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin));
811     fN2 = 0;
812     return kFALSE;
813   }
814
815   // update used clusters and select
816   fNUsed = 0;
817   for (Int_t it = 0; it < AliTRDtrackerV1::GetNTimeBins(); it++) {
818     if(fClusters[it] && fClusters[it]->IsUsed()) fNUsed++;
819     if(fClusters[it+35] && fClusters[it+35]->IsUsed()) fNUsed++;
820   }
821   if (fN2-fNUsed < kClmin){
822     //AliWarning(Form("Too many clusters already in use %d (from %d).", fNUsed, fN2));
823     fN2 = 0;
824     return kFALSE;
825   }
826
827   // set the Lorentz angle for this tracklet  
828   SetExB();
829
830   // calculate dx for time bins in the drift region (calibration aware)
831   Int_t irp = 0; Float_t x[2]={0., 0.}; Int_t tb[2] = {0, 0};
832   for (Int_t it = t0; it < AliTRDtrackerV1::GetNTimeBins(); it++) {
833     if(!fClusters[it]) continue;
834     x[irp]  = fClusters[it]->GetX();
835     tb[irp] = it;
836     irp++;
837     if(irp==2) break;
838   }
839   Int_t dtb = tb[1] - tb[0];
840   fdX = dtb ? (x[0] - x[1]) / dtb : 0.15;
841
842   // update X0 from the clusters (calibration/alignment aware) TODO remove dependence on x0 !!
843   for (Int_t it = 0; it < AliTRDtrackerV1::GetNTimeBins(); it++) {
844     if(!(layer = chamber->GetTB(it))) continue;
845     if(!layer->IsT0()) continue;
846     if(fClusters[it]){ 
847       fX0 = fClusters[it]->GetX();
848       break;
849     } else { // we have to infere the position of the anode wire from the other clusters
850       for (Int_t jt = it+1; jt < AliTRDtrackerV1::GetNTimeBins(); jt++) {
851         if(!fClusters[jt]) continue;
852         fX0 = fClusters[jt]->GetX() + fdX * (jt - it);
853         break;
854       }
855     }
856   }     
857
858   return kTRUE;
859 }
860
861 //____________________________________________________________
862 void AliTRDseedV1::Bootstrap(const AliTRDReconstructor *rec)
863 {
864 //   Fill in all derived information. It has to be called after recovery from file or HLT.
865 //   The primitive data are
866 //   - list of clusters
867 //   - detector (as the detector will be removed from clusters)
868 //   - position of anode wire (fX0) - temporary
869 //   - track reference position and direction
870 //   - momentum of the track
871 //   - time bin length [cm]
872 // 
873 //   A.Bercuci <A.Bercuci@gsi.de> Oct 30th 2008
874 //
875   fReconstructor = rec;
876   AliTRDgeometry g;
877   AliTRDpadPlane *pp = g.GetPadPlane(fDet);
878   fTilt      = TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle());
879   fPadLength = pp->GetLengthIPad();
880   fSnp = fYref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);
881   fTgl = fZref[1];
882   fN = 0; fN2 = 0; fMPads = 0.;
883   AliTRDcluster **cit = &fClusters[0];
884   for(Int_t ic = knTimebins; ic--; cit++){
885     if(!(*cit)) return;
886     fN++; fN2++;
887     fX[ic] = (*cit)->GetX() - fX0;
888     fY[ic] = (*cit)->GetY();
889     fZ[ic] = (*cit)->GetZ();
890   }
891   Update(); // Fit();
892   CookLabels();
893   GetProbability();
894 }
895
896
897 //____________________________________________________________________
898 Bool_t AliTRDseedV1::Fit(Bool_t tilt, Int_t errors)
899 {
900   //
901   // Linear fit of the tracklet
902   //
903   // Parameters :
904   //
905   // Output :
906   //  True if successful
907   //
908   // Detailed description
909   // 2. Check if tracklet crosses pad row boundary
910   // 1. Calculate residuals in the y (r-phi) direction
911   // 3. Do a Least Square Fit to the data
912   //
913
914   const Int_t kClmin = 8;
915
916
917   // cluster error parametrization parameters 
918   // 1. sy total charge
919   const Float_t sq0inv = 0.019962; // [1/q0]
920   const Float_t sqb    = 1.0281564;    //[cm]
921   // 2. sy for the PRF
922   const Float_t scy[AliTRDgeometry::kNlayer][4] = {
923     {2.827e-02, 9.600e-04, 4.296e-01, 2.271e-02},
924     {2.952e-02,-2.198e-04, 4.146e-01, 2.339e-02},
925     {3.090e-02, 1.514e-03, 4.020e-01, 2.402e-02},
926     {3.260e-02,-2.037e-03, 3.946e-01, 2.509e-02},
927     {3.439e-02,-3.601e-04, 3.883e-01, 2.623e-02},
928     {3.510e-02, 2.066e-03, 3.651e-01, 2.588e-02},
929   };
930   // 3. sy parallel to the track
931   const Float_t sy0 =  2.649e-02; // [cm]
932   const Float_t sya = -8.864e-04; // [cm]
933   const Float_t syb = -2.435e-01; // [cm]
934
935   // 4. sx parallel to the track
936   const Float_t sxgc = 5.427e-02;
937   const Float_t sxgm = 7.783e-01;
938   const Float_t sxgs = 2.743e-01;
939   const Float_t sxe0 =-2.065e+00;
940   const Float_t sxe1 =-2.978e-02;
941
942   // 5. sx perpendicular to the track
943 //   const Float_t sxd0 = 1.881e-02;
944 //   const Float_t sxd1 =-4.101e-01;
945 //   const Float_t sxd2 = 1.572e+00;
946
947   // get track direction
948   Double_t y0   = fYref[0];
949   Double_t dydx = fYref[1]; 
950   Double_t z0   = fZref[0];
951   Double_t dzdx = fZref[1];
952   Double_t yt, zt;
953
954   const Int_t kNtb = AliTRDtrackerV1::GetNTimeBins();
955   //AliTRDtrackerV1::AliTRDLeastSquare fitterZ;
956   TLinearFitter  fitterY(1, "pol1");
957   // convertion factor from square to gauss distribution for sigma
958   //Double_t convert = 1./TMath::Sqrt(12.);
959   
960   // book cluster information
961   Double_t q, xc[knTimebins], yc[knTimebins], zc[knTimebins], sy[knTimebins]/*, sz[knTimebins]*/;
962 //   Int_t zRow[knTimebins];
963   
964   Int_t ily = AliTRDgeometry::GetLayer(fDet);
965   fN = 0; //fXref = 0.; Double_t ssx = 0.;
966   AliTRDcluster *c=0x0, **jc = &fClusters[0];
967   for (Int_t ic=0; ic<kNtb; ic++, ++jc) {
968     //zRow[ic] = -1;
969     xc[ic]  = -1.;
970     yc[ic]  = 999.;
971     zc[ic]  = 999.;
972     sy[ic]  = 0.;
973     //sz[ic]  = 0.;
974     if(!(c = (*jc))) continue;
975     if(!c->IsInChamber()) continue;
976
977     Float_t w = 1.;
978     if(c->GetNPads()>4) w = .5;
979     if(c->GetNPads()>5) w = .2;
980
981     //zRow[fN] = c->GetPadRow();
982     // correct cluster position for PRF and v drift
983     Double_t x, y; GetClusterXY(c, x, y);
984     xc[fN]   = fX0 - x;
985     yc[fN]   = y;
986     zc[fN]   = c->GetZ();
987
988     // extrapolated y value for the track
989     yt = y0 - xc[fN]*dydx; 
990     // extrapolated z value for the track
991     zt = z0 - xc[fN]*dzdx; 
992     // tilt correction
993     if(tilt) yc[fN] -= fTilt*(zc[fN] - zt); 
994
995     // ELABORATE CLUSTER ERROR
996     // TODO to be moved to AliTRDcluster
997     q = TMath::Abs(c->GetQ());
998     Double_t tgg = (dydx-fExB)/(1.+dydx*fExB); tgg *= tgg;
999     // basic y error (|| to track).
1000     sy[fN]  = xc[fN] < AliTRDgeometry::CamHght() ? 2. : sy0 + sya*TMath::Exp(1./(xc[fN]+syb));
1001     //printf("cluster[%d]\n\tsy[0] = %5.3e [um]\n", fN,  sy[fN]*1.e4);
1002     // y error due to total charge
1003     sy[fN] += sqb*(1./q - sq0inv);
1004     //printf("\tsy[1] = %5.3e [um]\n", sy[fN]*1.e4);
1005     // y error due to PRF
1006     sy[fN] += scy[ily][0]*TMath::Gaus(c->GetCenter(), scy[ily][1], scy[ily][2]) - scy[ily][3];
1007     //printf("\tsy[2] = %5.3e [um]\n", sy[fN]*1.e4);
1008
1009     sy[fN] *= sy[fN];
1010
1011     // ADD ERROR ON x
1012     // error of drift length parallel to the track
1013     Double_t sx = sxgc*TMath::Gaus(xc[fN], sxgm, sxgs) + TMath::Exp(sxe0+sxe1*xc[fN]); // [cm]
1014     //printf("\tsx[0] = %5.3e [um]\n", sx*1.e4);
1015     // error of drift length perpendicular to the track
1016     //sx += sxd0 + sxd1*d + sxd2*d*d;
1017     sx *= sx; // square sx
1018     // update xref
1019     //fXref += xc[fN]/sx; ssx+=1./sx;
1020
1021     // add error from ExB 
1022     if(errors>0) sy[fN] += fExB*fExB*sx;
1023     //printf("\tsy[3] = %5.3e [um^2]\n", sy[fN]*1.e8);
1024
1025     // global radial error due to misalignment/miscalibration
1026     Double_t sx0  = 0.; sx0 *= sx0;
1027     // add sx contribution to sy due to track angle
1028     if(errors>1) sy[fN] += tgg*(sx+sx0);
1029     // TODO we should add tilt pad correction here
1030     //printf("\tsy[4] = %5.3e [um^2]\n", sy[fN]*1.e8);
1031     c->SetSigmaY2(sy[fN]);
1032
1033     sy[fN]  = TMath::Sqrt(sy[fN]);
1034     fitterY.AddPoint(&xc[fN], yc[fN]/*-yt*/, sy[fN]);
1035     fN++;
1036   }
1037   // to few clusters
1038   if (fN < kClmin) return kFALSE; 
1039
1040   // fit XY
1041   fitterY.Eval();
1042   fYfit[0] = fitterY.GetParameter(0);
1043   fYfit[1] = -fitterY.GetParameter(1);
1044   // store covariance
1045   Double_t *p = fitterY.GetCovarianceMatrix();
1046   fCov[0] = p[0]; // variance of y0
1047   fCov[1] = p[1]; // covariance of y0, dydx
1048   fCov[2] = p[3]; // variance of dydx
1049   // the ref radial position is set at the minimum of 
1050   // the y variance of the tracklet
1051   fXref = -fCov[1]/fCov[2]; //fXref = fX0 - fXref;
1052
1053   // fit XZ
1054   if(IsRowCross()){ 
1055     // TODO pad row cross position estimation !!!
1056     //AliInfo(Form("Padrow cross in detector %d", fDet));
1057     fZfit[0] = .5*(zc[0]+zc[fN-1]); fZfit[1] = 0.;
1058   } else {
1059     fZfit[0] = zc[0]; fZfit[1] = 0.;
1060   }
1061
1062
1063 //   // determine z offset of the fit
1064 //   Float_t zslope = 0.;
1065 //   Int_t nchanges = 0, nCross = 0;
1066 //   if(nz==2){ // tracklet is crossing pad row
1067 //     // Find the break time allowing one chage on pad-rows
1068 //     // with maximal number of accepted clusters
1069 //     Int_t padRef = zRow[0];
1070 //     for (Int_t ic=1; ic<fN; ic++) {
1071 //       if(zRow[ic] == padRef) continue;
1072 //       
1073 //       // debug
1074 //       if(zRow[ic-1] == zRow[ic]){
1075 //         printf("ERROR in pad row change!!!\n");
1076 //       }
1077 //     
1078 //       // evaluate parameters of the crossing point
1079 //       Float_t sx = (xc[ic-1] - xc[ic])*convert;
1080 //       fCross[0] = .5 * (xc[ic-1] + xc[ic]);
1081 //       fCross[2] = .5 * (zc[ic-1] + zc[ic]);
1082 //       fCross[3] = TMath::Max(dzdx * sx, .01);
1083 //       zslope    = zc[ic-1] > zc[ic] ? 1. : -1.;
1084 //       padRef    = zRow[ic];
1085 //       nCross    = ic;
1086 //       nchanges++;
1087 //     }
1088 //   }
1089 // 
1090 //   // condition on nCross and reset nchanges TODO
1091 // 
1092 //   if(nchanges==1){
1093 //     if(dzdx * zslope < 0.){
1094 //       AliInfo("Tracklet-Track mismatch in dzdx. TODO.");
1095 //     }
1096 // 
1097 // 
1098 //     //zc[nc] = fitterZ.GetFunctionParameter(0); 
1099 //     fCross[1] = fYfit[0] - fCross[0] * fYfit[1];
1100 //     fCross[0] = fX0 - fCross[0];
1101 //   }
1102
1103   UpdateUsed();
1104   return kTRUE;
1105 }
1106
1107
1108 //___________________________________________________________________
1109 void AliTRDseedV1::Print(Option_t *o) const
1110 {
1111   //
1112   // Printing the seedstatus
1113   //
1114
1115   AliInfo(Form("Det[%3d] Tilt[%+6.2f] Pad[%5.2f]", fDet, fTilt, fPadLength));
1116   AliInfo(Form("Nattach[%2d] Nfit[%2d] Nuse[%2d] pads[%f]", fN, fN2, fNUsed, fMPads));
1117   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]));
1118   AliInfo(Form("Ref        y[%7.2f] z[%7.2f] dydx[%5.2f] dzdx[%5.2f]", fYref[0], fZref[0], fYref[1], fZref[1]))
1119
1120
1121   if(strcmp(o, "a")!=0) return;
1122
1123   AliTRDcluster* const* jc = &fClusters[0];
1124   for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++, jc++) {
1125     if(!(*jc)) continue;
1126     (*jc)->Print(o);
1127   }
1128 }
1129
1130
1131 //___________________________________________________________________
1132 Bool_t AliTRDseedV1::IsEqual(const TObject *o) const
1133 {
1134   // Checks if current instance of the class has the same essential members
1135   // as the given one
1136
1137   if(!o) return kFALSE;
1138   const AliTRDseedV1 *inTracklet = dynamic_cast<const AliTRDseedV1*>(o);
1139   if(!inTracklet) return kFALSE;
1140
1141   for (Int_t i = 0; i < 2; i++){
1142     if ( fYref[i] != inTracklet->GetYref(i) ) return kFALSE;
1143     if ( fZref[i] != inTracklet->GetZref(i) ) return kFALSE;
1144   }
1145   
1146   if ( fSigmaY != inTracklet->GetSigmaY() ) return kFALSE;
1147   if ( fSigmaY2 != inTracklet->GetSigmaY2() ) return kFALSE;
1148   if ( fTilt != inTracklet->GetTilt() ) return kFALSE;
1149   if ( fPadLength != inTracklet->GetPadLength() ) return kFALSE;
1150   
1151   for (Int_t i = 0; i < knTimebins; i++){
1152     if ( fX[i] != inTracklet->GetX(i) ) return kFALSE;
1153     if ( fY[i] != inTracklet->GetY(i) ) return kFALSE;
1154     if ( fZ[i] != inTracklet->GetZ(i) ) return kFALSE;
1155     if ( fIndexes[i] != inTracklet->GetIndexes(i) ) return kFALSE;
1156     if ( fUsable[i] != inTracklet->IsUsable(i) ) return kFALSE;
1157   }
1158
1159   for (Int_t i=0; i < 2; i++){
1160     if ( fYfit[i] != inTracklet->GetYfit(i) ) return kFALSE;
1161     if ( fZfit[i] != inTracklet->GetZfit(i) ) return kFALSE;
1162     if ( fYfitR[i] != inTracklet->GetYfitR(i) ) return kFALSE;
1163     if ( fZfitR[i] != inTracklet->GetZfitR(i) ) return kFALSE;
1164     if ( fLabels[i] != inTracklet->GetLabels(i) ) return kFALSE;
1165   }
1166   
1167   if ( fMeanz != inTracklet->GetMeanz() ) return kFALSE;
1168   if ( fZProb != inTracklet->GetZProb() ) return kFALSE;
1169   if ( fN2 != inTracklet->GetN2() ) return kFALSE;
1170   if ( fNUsed != inTracklet->GetNUsed() ) return kFALSE;
1171   if ( fFreq != inTracklet->GetFreq() ) return kFALSE;
1172   if ( fNChange != inTracklet->GetNChange() ) return kFALSE;
1173   if ( fNChange != inTracklet->GetNChange() ) return kFALSE;
1174    
1175   if ( fC != inTracklet->GetC() ) return kFALSE;
1176   if ( fCC != inTracklet->GetCC() ) return kFALSE;
1177   if ( fChi2 != inTracklet->GetChi2() ) return kFALSE;
1178   //  if ( fChi2Z != inTracklet->GetChi2Z() ) return kFALSE;
1179
1180   if ( fDet != inTracklet->GetDetector() ) return kFALSE;
1181   if ( fMom != inTracklet->GetMomentum() ) return kFALSE;
1182   if ( fdX != inTracklet->GetdX() ) return kFALSE;
1183   
1184   for (Int_t iCluster = 0; iCluster < knTimebins; iCluster++){
1185     AliTRDcluster *curCluster = fClusters[iCluster];
1186     AliTRDcluster *inCluster = inTracklet->GetClusters(iCluster);
1187     if (curCluster && inCluster){
1188       if (! curCluster->IsEqual(inCluster) ) {
1189         curCluster->Print();
1190         inCluster->Print();
1191         return kFALSE;
1192       }
1193     } else {
1194       // if one cluster exists, and corresponding 
1195       // in other tracklet doesn't - return kFALSE
1196       if(curCluster || inCluster) return kFALSE;
1197     }
1198   }
1199   return kTRUE;
1200 }