added track covariance matrix (yz plane) to the tracklet
[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     (*AliTRDtrackerV1::DebugStreamer()) << "AttachClustersIter"
388       << "chamber.="   << &ch
389       << "tracklet.="  << this
390       << "\n";  
391   }
392
393   Float_t  tquality;
394   Double_t kroady = fReconstructor->GetRecoParam() ->GetRoad1y();
395   Double_t kroadz = fPadLength * .5 + 1.;
396   
397   // initialize configuration parameters
398   Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
399   Int_t   niter = kZcorr ? 1 : 2;
400   
401   Double_t yexp, zexp;
402   Int_t ncl = 0;
403   // start seed update
404   for (Int_t iter = 0; iter < niter; iter++) {
405     ncl = 0;
406     for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
407       if(!(layer = chamber->GetTB(iTime))) continue;
408       if(!Int_t(*layer)) continue;
409       
410       // define searching configuration
411       Double_t dxlayer = layer->GetX() - fX0;
412       if(c){
413         zexp = c->GetZ();
414         //Try 2 pad-rows in second iteration
415         if (iter > 0) {
416           zexp = fZref[0] + fZref[1] * dxlayer - zcorr;
417           if (zexp > c->GetZ()) zexp = c->GetZ() + fPadLength*0.5;
418           if (zexp < c->GetZ()) zexp = c->GetZ() - fPadLength*0.5;
419         }
420       } else zexp = fZref[0] + (kZcorr ? fZref[1] * dxlayer : 0.);
421       yexp  = fYref[0] + fYref[1] * dxlayer - zcorr;
422       
423       // Get and register cluster
424       Int_t    index = layer->SearchNearestCluster(yexp, zexp, kroady, kroadz);
425       if (index < 0) continue;
426       AliTRDcluster *cl = (*layer)[index];
427       
428       fIndexes[iTime]  = layer->GetGlobalIndex(index);
429       fClusters[iTime] = cl;
430       fY[iTime]        = cl->GetY();
431       fZ[iTime]        = cl->GetZ();
432       ncl++;
433     }
434     if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d ncl [%d] = %d", iter, fDet, ncl));
435     
436     if(ncl>1){  
437       // calculate length of the time bin (calibration aware)
438       Int_t irp = 0; Float_t x[2]; Int_t tb[2];
439       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
440         if(!fClusters[iTime]) continue;
441         x[irp]  = fClusters[iTime]->GetX();
442         tb[irp] = iTime;
443         irp++;
444         if(irp==2) break;
445       } 
446       fdX = (x[1] - x[0]) / (tb[0] - tb[1]);
447   
448       // update X0 from the clusters (calibration/alignment aware)
449       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
450         if(!(layer = chamber->GetTB(iTime))) continue;
451         if(!layer->IsT0()) continue;
452         if(fClusters[iTime]){ 
453           fX0 = fClusters[iTime]->GetX();
454           break;
455         } else { // we have to infere the position of the anode wire from the other clusters
456           for (Int_t jTime = iTime+1; jTime < AliTRDtrackerV1::GetNTimeBins(); jTime++) {
457             if(!fClusters[jTime]) continue;
458             fX0 = fClusters[jTime]->GetX() + fdX * (jTime - iTime);
459             break;
460           }
461         }
462       } 
463       
464       // update YZ reference point
465       // TODO
466       
467       // update x reference positions (calibration/alignment aware)
468       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
469         if(!fClusters[iTime]) continue;
470         fX[iTime] = fClusters[iTime]->GetX() - fX0;
471       } 
472       
473       AliTRDseed::Update();
474     }
475     if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d nclFit [%d] = %d", iter, fDet, fN2));
476     
477     if(IsOK()){
478       tquality = GetQuality(kZcorr);
479       if(tquality < quality) break;
480       else quality = tquality;
481     }
482     kroadz *= 2.;
483   } // Loop: iter
484   if (!IsOK()) return kFALSE;
485
486   if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=1) CookLabels();
487   UpdateUsed();
488   return kTRUE; 
489 }
490
491 //____________________________________________________________________
492 Bool_t  AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *chamber
493                                       ,Bool_t kZcorr)
494 {
495   //
496   // Projective algorithm to attach clusters to seeding tracklets
497   //
498   // Parameters
499   //
500   // Output
501   //
502   // Detailed description
503   // 1. Collapse x coordinate for the full detector plane
504   // 2. truncated mean on y (r-phi) direction
505   // 3. purge clusters
506   // 4. truncated mean on z direction
507   // 5. purge clusters
508   // 6. fit tracklet
509   //    
510
511   if(!fReconstructor->GetRecoParam() ){
512     AliError("Seed can not be used without a valid RecoParam.");
513     return kFALSE;
514   }
515
516   const Int_t kClusterCandidates = 2 * knTimebins;
517   
518   //define roads
519   Double_t kroady = fReconstructor->GetRecoParam() ->GetRoad1y();
520   Double_t kroadz = fPadLength * 1.5 + 1.;
521   // correction to y for the tilting angle
522   Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
523
524   // working variables
525   AliTRDcluster *clusters[kClusterCandidates];
526   Double_t cond[4], yexp[knTimebins], zexp[knTimebins],
527     yres[kClusterCandidates], zres[kClusterCandidates];
528   Int_t ncl, *index = 0x0, tboundary[knTimebins];
529   
530   // Do cluster projection
531   AliTRDchamberTimeBin *layer = 0x0;
532   Int_t nYclusters = 0; Bool_t kEXIT = kFALSE;
533   for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
534     if(!(layer = chamber->GetTB(iTime))) continue;
535     if(!Int_t(*layer)) continue;
536     
537     fX[iTime] = layer->GetX() - fX0;
538     zexp[iTime] = fZref[0] + fZref[1] * fX[iTime];
539     yexp[iTime] = fYref[0] + fYref[1] * fX[iTime] - zcorr;
540     
541     // build condition and process clusters
542     cond[0] = yexp[iTime] - kroady; cond[1] = yexp[iTime] + kroady;
543     cond[2] = zexp[iTime] - kroadz; cond[3] = zexp[iTime] + kroadz;
544     layer->GetClusters(cond, index, ncl);
545     for(Int_t ic = 0; ic<ncl; ic++){
546       AliTRDcluster *c = layer->GetCluster(index[ic]);
547       clusters[nYclusters] = c;
548       yres[nYclusters++] = c->GetY() - yexp[iTime];
549       if(nYclusters >= kClusterCandidates) {
550         AliWarning(Form("Cluster candidates reached limit %d. Some may be lost.", kClusterCandidates));
551         kEXIT = kTRUE;
552         break;
553       }
554     }
555     tboundary[iTime] = nYclusters;
556     if(kEXIT) break;
557   }
558   
559   // Evaluate truncated mean on the y direction
560   Double_t mean, sigma;
561   AliMathBase::EvaluateUni(nYclusters, yres, mean, sigma, Int_t(nYclusters*.8)-2);
562   // purge cluster candidates
563   Int_t nZclusters = 0;
564   for(Int_t ic = 0; ic<nYclusters; ic++){
565     if(yres[ic] - mean > 4. * sigma){
566       clusters[ic] = 0x0;
567       continue;
568     }
569     zres[nZclusters++] = clusters[ic]->GetZ() - zexp[clusters[ic]->GetLocalTimeBin()];
570   }
571   
572   // Evaluate truncated mean on the z direction
573   AliMathBase::EvaluateUni(nZclusters, zres, mean, sigma, Int_t(nZclusters*.8)-2);
574   // purge cluster candidates
575   for(Int_t ic = 0; ic<nZclusters; ic++){
576     if(zres[ic] - mean > 4. * sigma){
577       clusters[ic] = 0x0;
578       continue;
579     }
580   }
581
582   
583   // Select only one cluster/TimeBin
584   Int_t lastCluster = 0;
585   fN2 = 0;
586   for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
587     ncl = tboundary[iTime] - lastCluster;
588     if(!ncl) continue;
589     Int_t iptr = lastCluster;
590     if(ncl > 1){
591       Float_t dold = 9999.;
592       for(int ic=lastCluster; ic<tboundary[iTime]; ic++){
593         if(!clusters[ic]) continue;
594         Float_t y = yexp[iTime] - clusters[ic]->GetY();
595         Float_t z = zexp[iTime] - clusters[ic]->GetZ();
596         Float_t d = y * y + z * z;
597         if(d > dold) continue;
598         dold = d;
599         iptr = ic;
600       }
601     }
602     fIndexes[iTime]  = chamber->GetTB(iTime)->GetGlobalIndex(iptr);
603     fClusters[iTime] = clusters[iptr];
604     fY[iTime]        = clusters[iptr]->GetY();
605     fZ[iTime]        = clusters[iptr]->GetZ();
606     lastCluster      = tboundary[iTime];
607     fN2++;
608   }
609   
610   // number of minimum numbers of clusters expected for the tracklet
611   Int_t kClmin = Int_t(fReconstructor->GetRecoParam() ->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
612   if (fN2 < kClmin){
613     AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin));
614     fN2 = 0;
615     return kFALSE;
616   }
617
618   // update used clusters
619   fNUsed = 0;
620   for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
621     if(!fClusters[iTime]) continue;
622     if((fClusters[iTime]->IsUsed())) fNUsed++;
623   }
624
625   if (fN2-fNUsed < kClmin){
626     AliWarning(Form("Too many clusters already in use %d (from %d).", fNUsed, fN2));
627     fN2 = 0;
628     return kFALSE;
629   }
630   
631   return kTRUE;
632 }
633
634 //____________________________________________________________
635 void AliTRDseedV1::Bootstrap(const AliTRDReconstructor *rec)
636 {
637 //   Fill in all derived information. It has to be called after recovery from file or HLT.
638 //   The primitive data are
639 //   - list of clusters
640 //   - detector (as the detector will be removed from clusters)
641 //   - position of anode wire (fX0) - temporary
642 //   - track reference position and direction
643 //   - momentum of the track
644 //   - time bin length [cm]
645 // 
646 //   A.Bercuci <A.Bercuci@gsi.de> Oct 30th 2008
647 //
648   fReconstructor = rec;
649   AliTRDgeometry g;
650   AliTRDpadPlane *pp = g.GetPadPlane(fDet);
651   fTilt      = TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle());
652   fPadLength = pp->GetLengthIPad();
653   fSnp = fYref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);
654   fTgl = fZref[1];
655   fN = 0; fN2 = 0; fMPads = 0.;
656   AliTRDcluster **cit = &fClusters[0];
657   for(Int_t ic = knTimebins; ic--; cit++){
658     if(!(*cit)) return;
659     fN++; fN2++;
660     fX[ic] = (*cit)->GetX() - fX0;
661     fY[ic] = (*cit)->GetY();
662     fZ[ic] = (*cit)->GetZ();
663   }
664   Update(); // Fit();
665   CookLabels();
666   GetProbability();
667 }
668
669
670 //____________________________________________________________________
671 Bool_t AliTRDseedV1::Fit(Bool_t tilt)
672 {
673   //
674   // Linear fit of the tracklet
675   //
676   // Parameters :
677   //
678   // Output :
679   //  True if successful
680   //
681   // Detailed description
682   // 2. Check if tracklet crosses pad row boundary
683   // 1. Calculate residuals in the y (r-phi) direction
684   // 3. Do a Least Square Fit to the data
685   //
686
687   const Int_t kClmin = 8;
688   const Float_t q0 = 100.;
689   const Float_t clSigma0 = 2.E-2;    //[cm]
690   const Float_t clSlopeQ = -1.19E-2; //[1/cm]
691
692   // get track direction
693   Double_t y0   = fYref[0];
694   Double_t dydx = fYref[1]; 
695   Double_t z0   = fZref[0];
696   Double_t dzdx = fZref[1];
697   Double_t yt, zt;
698
699   const Int_t kNtb = AliTRDtrackerV1::GetNTimeBins();
700   AliTRDtrackerV1::AliTRDLeastSquare fitterY, fitterZ;
701
702   // convertion factor from square to gauss distribution for sigma
703   Double_t convert = 1./TMath::Sqrt(12.);
704   
705   // book cluster information
706   Double_t xc[knTimebins], yc[knTimebins], zc[knTimebins], sy[knTimebins], sz[knTimebins];
707   Int_t zRow[knTimebins];
708
709
710   fN = 0;
711   AliTRDcluster *c=0x0, **jc = &fClusters[0];
712   for (Int_t ic=0; ic<kNtb; ic++, ++jc) {
713     zRow[ic] = -1;
714     xc[ic]  = -1.;
715     yc[ic]  = 999.;
716     zc[ic]  = 999.;
717     sy[ic]  = 0.;
718     sz[ic]  = 0.;
719     if(!(c = (*jc))) continue;
720     if(!c->IsInChamber()) continue;
721     Float_t w = 1.;
722     if(c->GetNPads()>4) w = .5;
723     if(c->GetNPads()>5) w = .2;
724     zRow[fN] = c->GetPadRow();
725     xc[fN]   = fX0 - c->GetX();
726     yc[fN]   = c->GetY();
727     zc[fN]   = c->GetZ();
728
729     // extrapolated y value for the track
730     yt = y0 - xc[fN]*dydx; 
731     // extrapolated z value for the track
732     zt = z0 - xc[fN]*dzdx; 
733     // tilt correction
734     if(tilt) yc[fN] -= fTilt*(zc[fN] - zt); 
735
736     // elaborate cluster error
737     Float_t qr = c->GetQ() - q0;
738     sy[fN]   = qr < 0. ? clSigma0*TMath::Exp(clSlopeQ*qr) : clSigma0;
739
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.GetFunctionParameter(0);
752   fYfit[1] = dydx-fitterY.GetFunctionParameter(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 }