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