Major update for the TRD tracking code
[u/mrichter/AliRoot.git] / TRD / AliTRDseedV1.cxx
1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 *                                                                        *
4 * Author: The ALICE Off-line Project.                                    *
5 * Contributors are mentioned in the code where appropriate.              *
6 *                                                                        *
7 * Permission to use, copy, modify and distribute this software and its   *
8 * documentation strictly for non-commercial purposes is hereby granted   *
9 * without fee, provided that the above copyright notice appears in all   *
10 * copies and that both the copyright notice and this permission notice   *
11 * appear in the supporting documentation. The authors make no claims     *
12 * about the suitability of this software for any purpose. It is          *
13 * provided "as is" without express or implied warranty.                  *
14 **************************************************************************/
15
16 /* $Id$ */
17
18 ////////////////////////////////////////////////////////////////////////////
19 //                                                                        //
20 //  The TRD track seed                                                    //
21 //                                                                        //
22 //  Authors:                                                              //
23 //    Alex Bercuci <A.Bercuci@gsi.de>                                     //
24 //    Markus Fasel <M.Fasel@gsi.de>                                       //
25 //                                                                        //
26 ////////////////////////////////////////////////////////////////////////////
27
28 #include "TMath.h"
29 #include "TLinearFitter.h"
30 #include "TClonesArray.h" // tmp
31 #include <TTreeStream.h>
32
33 #include "AliLog.h"
34 #include "AliMathBase.h"
35 #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   :TObject()
59   ,fReconstructor(0x0)
60   ,fClusterIter(0x0)
61   ,fExB(0.)
62   ,fVD(0.)
63   ,fT0(0.)
64   ,fS2PRF(0.)
65   ,fDiffL(0.)
66   ,fDiffT(0.)
67   ,fClusterIdx(0)
68   ,fUsable(0)
69   ,fN2(0)
70   ,fNUsed(0)
71   ,fDet(det)
72   ,fTilt(0.)
73   ,fPadLength(0.)
74   ,fMom(0.)
75   ,fdX(0.)
76   ,fX0(0.)
77   ,fX(0.)
78   ,fY(0.)
79   ,fZ(0.)
80   ,fS2Y(0.)
81   ,fS2Z(0.)
82   ,fC(0.)
83   ,fChi2(0.)
84 {
85   //
86   // Constructor
87   //
88   for(Int_t ic=kNTimeBins; ic--;) fIndexes[ic] = -1;
89   memset(fClusters, 0, kNTimeBins*sizeof(AliTRDcluster*));
90   fYref[0] = 0.; fYref[1] = 0.; 
91   fZref[0] = 0.; fZref[1] = 0.; 
92   fYfit[0] = 0.; fYfit[1] = 0.; 
93   fZfit[0] = 0.; fZfit[1] = 0.; 
94   memset(fdEdx, 0, kNSlices*sizeof(Float_t)); 
95   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec]  = -1.;
96   fLabels[0]=-1; fLabels[1]=-1; // most freq MC labels
97   fLabels[2]=0;  // number of different labels for tracklet
98   fRefCov[0] = 1.; fRefCov[1] = 0.; fRefCov[2] = 1.;
99   // covariance matrix [diagonal]
100   // default sy = 200um and sz = 2.3 cm 
101   fCov[0] = 4.e-4; fCov[1] = 0.; fCov[2] = 5.3; 
102 }
103
104 //____________________________________________________________________
105 AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
106   :TObject((TObject&)ref)
107   ,fReconstructor(0x0)
108   ,fClusterIter(0x0)
109   ,fExB(0.)
110   ,fVD(0.)
111   ,fT0(0.)
112   ,fS2PRF(0.)
113   ,fDiffL(0.)
114   ,fDiffT(0.)
115   ,fClusterIdx(0)
116   ,fUsable(0)
117   ,fN2(0)
118   ,fNUsed(0)
119   ,fDet(-1)
120   ,fTilt(0.)
121   ,fPadLength(0.)
122   ,fMom(0.)
123   ,fdX(0.)
124   ,fX0(0.)
125   ,fX(0.)
126   ,fY(0.)
127   ,fZ(0.)
128   ,fS2Y(0.)
129   ,fS2Z(0.)
130   ,fC(0.)
131   ,fChi2(0.)
132 {
133   //
134   // Copy Constructor performing a deep copy
135   //
136   if(this != &ref){
137     ref.Copy(*this);
138   }
139   SetBit(kOwner, kFALSE);
140 }
141
142
143 //____________________________________________________________________
144 AliTRDseedV1& AliTRDseedV1::operator=(const AliTRDseedV1 &ref)
145 {
146   //
147   // Assignment Operator using the copy function
148   //
149
150   if(this != &ref){
151     ref.Copy(*this);
152   }
153   SetBit(kOwner, kFALSE);
154
155   return *this;
156 }
157
158 //____________________________________________________________________
159 AliTRDseedV1::~AliTRDseedV1()
160 {
161   //
162   // Destructor. The RecoParam object belongs to the underlying tracker.
163   //
164
165   //printf("I-AliTRDseedV1::~AliTRDseedV1() : Owner[%s]\n", IsOwner()?"YES":"NO");
166
167   if(IsOwner()) {
168     for(int itb=0; itb<kNTimeBins; itb++){
169       if(!fClusters[itb]) continue; 
170       //AliInfo(Form("deleting c %p @ %d", fClusters[itb], itb));
171       delete fClusters[itb];
172       fClusters[itb] = 0x0;
173     }
174   }
175 }
176
177 //____________________________________________________________________
178 void AliTRDseedV1::Copy(TObject &ref) const
179 {
180   //
181   // Copy function
182   //
183
184   //AliInfo("");
185   AliTRDseedV1 &target = (AliTRDseedV1 &)ref; 
186
187   target.fReconstructor = fReconstructor;
188   target.fClusterIter   = 0x0;
189   target.fExB           = fExB;
190   target.fVD            = fVD;
191   target.fT0            = fT0;
192   target.fS2PRF         = fS2PRF;
193   target.fDiffL         = fDiffL;
194   target.fDiffT         = fDiffT;
195   target.fClusterIdx    = 0;
196   target.fUsable        = fUsable;
197   target.fN2            = fN2;
198   target.fNUsed         = fNUsed;
199   target.fDet           = fDet;
200   target.fTilt          = fTilt;
201   target.fPadLength     = fPadLength;
202   target.fMom           = fMom;
203   target.fdX            = fdX;
204   target.fX0            = fX0;
205   target.fX             = fX;
206   target.fY             = fY;
207   target.fZ             = fZ;
208   target.fS2Y           = fS2Y;
209   target.fS2Z           = fS2Z;
210   target.fC             = fC;
211   target.fChi2          = fChi2;
212   
213   memcpy(target.fIndexes, fIndexes, kNTimeBins*sizeof(Int_t));
214   memcpy(target.fClusters, fClusters, kNTimeBins*sizeof(AliTRDcluster*));
215   target.fYref[0] = fYref[0]; target.fYref[1] = fYref[1]; 
216   target.fZref[0] = fZref[0]; target.fZref[1] = fZref[1]; 
217   target.fYfit[0] = fYfit[0]; target.fYfit[1] = fYfit[1]; 
218   target.fZfit[0] = fZfit[0]; target.fZfit[1] = fZfit[1]; 
219   memcpy(target.fdEdx, fdEdx, kNSlices*sizeof(Float_t)); 
220   memcpy(target.fProb, fProb, AliPID::kSPECIES*sizeof(Float_t)); 
221   memcpy(target.fLabels, fLabels, 3*sizeof(Int_t)); 
222   memcpy(target.fRefCov, fRefCov, 3*sizeof(Double_t)); 
223   memcpy(target.fCov, fCov, 3*sizeof(Double_t)); 
224   
225   TObject::Copy(ref);
226 }
227
228
229 //____________________________________________________________
230 Bool_t AliTRDseedV1::Init(AliTRDtrackV1 *track)
231 {
232 // Initialize this tracklet using the track information
233 //
234 // Parameters:
235 //   track - the TRD track used to initialize the tracklet
236 // 
237 // Detailed description
238 // The function sets the starting point and direction of the
239 // tracklet according to the information from the TRD track.
240 // 
241 // Caution
242 // The TRD track has to be propagated to the beginning of the
243 // chamber where the tracklet will be constructed
244 //
245
246   Double_t y, z; 
247   if(!track->GetProlongation(fX0, y, z)) return kFALSE;
248   UpDate(track);
249   return kTRUE;
250 }
251
252
253 //_____________________________________________________________________________
254 void AliTRDseedV1::Reset()
255 {
256   //
257   // Reset seed
258   //
259   fExB=0.;fVD=0.;fT0=0.;fS2PRF=0.;
260   fDiffL=0.;fDiffT=0.;
261   fClusterIdx=0;fUsable=0;
262   fN2=0;fNUsed=0;
263   fDet=-1;fTilt=0.;fPadLength=0.;
264   fMom=0.;
265   fdX=0.;fX0=0.; fX=0.; fY=0.; fZ=0.;
266   fS2Y=0.; fS2Z=0.;
267   fC=0.; fChi2 = 0.;
268
269   for(Int_t ic=kNTimeBins; ic--;) fIndexes[ic] = -1;
270   memset(fClusters, 0, kNTimeBins*sizeof(AliTRDcluster*));
271   fYref[0] = 0.; fYref[1] = 0.; 
272   fZref[0] = 0.; fZref[1] = 0.; 
273   fYfit[0] = 0.; fYfit[1] = 0.; 
274   fZfit[0] = 0.; fZfit[1] = 0.; 
275   memset(fdEdx, 0, kNSlices*sizeof(Float_t)); 
276   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec]  = -1.;
277   fLabels[0]=-1; fLabels[1]=-1; // most freq MC labels
278   fLabels[2]=0;  // number of different labels for tracklet
279   fRefCov[0] = 1.; fRefCov[1] = 0.; fRefCov[2] = 1.;
280   // covariance matrix [diagonal]
281   // default sy = 200um and sz = 2.3 cm 
282   fCov[0] = 4.e-4; fCov[1] = 0.; fCov[2] = 5.3; 
283 }
284
285 //____________________________________________________________________
286 void AliTRDseedV1::UpDate(const AliTRDtrackV1 *trk)
287
288   // update tracklet reference position from the TRD track
289   // Funny name to avoid the clash with the function AliTRDseed::Update() (has to be made obselete)
290
291   Double_t fSnp = trk->GetSnp();
292   Double_t fTgl = trk->GetTgl();
293   fMom = trk->GetP();
294   fYref[1] = fSnp/(1. - fSnp*fSnp);
295   fZref[1] = fTgl;
296   SetCovRef(trk->GetCovariance());
297
298   Double_t dx = trk->GetX() - fX0;
299   fYref[0] = trk->GetY() - dx*fYref[1];
300   fZref[0] = trk->GetZ() - dx*fZref[1];
301 }
302
303 //_____________________________________________________________________________
304 void AliTRDseedV1::UpdateUsed()
305 {
306   //
307   // Update used seed
308   //
309
310   fNUsed = 0;
311   for (Int_t i = kNTimeBins; i--; ) {
312     if (!fClusters[i]) continue;
313     if(!TESTBIT(fUsable, i)) continue;
314     if((fClusters[i]->IsUsed())) fNUsed++;
315   }
316 }
317
318 //_____________________________________________________________________________
319 void AliTRDseedV1::UseClusters()
320 {
321   //
322   // Use clusters
323   //
324   AliTRDcluster **c = &fClusters[0];
325   for (Int_t ic=kNTimeBins; ic--; c++) {
326     if(!(*c)) continue;
327     if(!((*c)->IsUsed())) (*c)->Use();
328   }
329 }
330
331
332 //____________________________________________________________________
333 void AliTRDseedV1::CookdEdx(Int_t nslices)
334 {
335 // Calculates average dE/dx for all slices and store them in the internal array fdEdx. 
336 //
337 // Parameters:
338 //  nslices : number of slices for which dE/dx should be calculated
339 // Output:
340 //  store results in the internal array fdEdx. This can be accessed with the method
341 //  AliTRDseedV1::GetdEdx()
342 //
343 // Detailed description
344 // Calculates average dE/dx for all slices. Depending on the PID methode 
345 // the number of slices can be 3 (LQ) or 8(NN). 
346 // The calculation of dQ/dl are done using the tracklet fit results (see AliTRDseedV1::GetdQdl(Int_t))
347 //
348 // The following effects are included in the calculation:
349 // 1. calibration values for t0 and vdrift (using x coordinate to calculate slice)
350 // 2. cluster sharing (optional see AliTRDrecoParam::SetClusterSharing())
351 // 3. cluster size
352 //
353
354   Int_t nclusters[kNSlices]; 
355   memset(nclusters, 0, kNSlices*sizeof(Int_t));
356   memset(fdEdx, 0, kNSlices*sizeof(Float_t));
357
358   const Double_t kDriftLength = (.5 * AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
359
360   AliTRDcluster *c = 0x0;
361   for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++){
362     if(!(c = fClusters[ic]) && !(c = fClusters[ic+kNtb])) continue;
363     Float_t dx = TMath::Abs(fX0 - c->GetX());
364     
365     // Filter clusters for dE/dx calculation
366     
367     // 1.consider calibration effects for slice determination
368     Int_t slice;
369     if(dx<kDriftLength){ // TODO should be replaced by c->IsInChamber() 
370       slice = Int_t(dx * nslices / kDriftLength);
371     } else slice = c->GetX() < fX0 ? nslices-1 : 0;
372
373
374     // 2. take sharing into account
375     Float_t w = c->IsShared() ? .5 : 1.;
376     
377     // 3. take into account large clusters TODO
378     //w *= c->GetNPads() > 3 ? .8 : 1.;
379     
380     //CHECK !!!
381     fdEdx[slice]   += w * GetdQdl(ic); //fdQdl[ic];
382     nclusters[slice]++;
383   } // End of loop over clusters
384
385   //if(fReconstructor->GetPIDMethod() == AliTRDReconstructor::kLQPID){
386   if(nslices == AliTRDpidUtil::kLQslices){
387   // calculate mean charge per slice (only LQ PID)
388     for(int is=0; is<nslices; is++){ 
389       if(nclusters[is]) fdEdx[is] /= nclusters[is];
390     }
391   }
392 }
393
394 //_____________________________________________________________________________
395 void AliTRDseedV1::CookLabels()
396 {
397   //
398   // Cook 2 labels for seed
399   //
400
401   Int_t labels[200];
402   Int_t out[200];
403   Int_t nlab = 0;
404   for (Int_t i = 0; i < kNTimeBins; i++) {
405     if (!fClusters[i]) continue;
406     for (Int_t ilab = 0; ilab < 3; ilab++) {
407       if (fClusters[i]->GetLabel(ilab) >= 0) {
408         labels[nlab] = fClusters[i]->GetLabel(ilab);
409         nlab++;
410       }
411     }
412   }
413
414   fLabels[2] = AliTRDtrackerV1::Freq(nlab,labels,out,kTRUE);
415   fLabels[0] = out[0];
416   if ((fLabels[2]  > 1) && (out[3] > 1)) fLabels[1] = out[2];
417 }
418
419
420 //____________________________________________________________________
421 void AliTRDseedV1::GetClusterXY(const AliTRDcluster *c, Double_t &x, Double_t &y)
422 {
423 // Return corrected position of the cluster taking into 
424 // account variation of the drift velocity with drift length.
425
426
427   // drift velocity correction TODO to be moved to the clusterizer
428   const Float_t cx[] = {
429     -9.6280e-02, 1.3091e-01,-1.7415e-02,-9.9221e-02,-1.2040e-01,-9.5493e-02,
430     -5.0041e-02,-1.6726e-02, 3.5756e-03, 1.8611e-02, 2.6378e-02, 3.3823e-02,
431      3.4811e-02, 3.5282e-02, 3.5386e-02, 3.6047e-02, 3.5201e-02, 3.4384e-02,
432      3.2864e-02, 3.1932e-02, 3.2051e-02, 2.2539e-02,-2.5154e-02,-1.2050e-01,
433     -1.2050e-01
434   };
435
436   // PRF correction TODO to be replaced by the gaussian 
437   // approximation with full error parametrization and // moved to the clusterizer
438   const Float_t cy[AliTRDgeometry::kNlayer][3] = {
439     { 4.014e-04, 8.605e-03, -6.880e+00},
440     {-3.061e-04, 9.663e-03, -6.789e+00},
441     { 1.124e-03, 1.105e-02, -6.825e+00},
442     {-1.527e-03, 1.231e-02, -6.777e+00},
443     { 2.150e-03, 1.387e-02, -6.783e+00},
444     {-1.296e-03, 1.486e-02, -6.825e+00}
445   }; 
446
447   Int_t ily = AliTRDgeometry::GetLayer(c->GetDetector());
448   x = c->GetX() - cx[c->GetLocalTimeBin()];
449   y = c->GetY() + cy[ily][0] + cy[ily][1] * TMath::Sin(cy[ily][2] * c->GetCenter());
450   return;
451 }
452
453 //____________________________________________________________________
454 Float_t AliTRDseedV1::GetdQdl(Int_t ic) const
455 {
456 // Using the linear approximation of the track inside one TRD chamber (TRD tracklet) 
457 // the charge per unit length can be written as:
458 // BEGIN_LATEX
459 // #frac{dq}{dl} = #frac{q_{c}}{dx * #sqrt{1 + #(){#frac{dy}{dx}}^{2}_{fit} + #(){#frac{dy}{dx}}^{2}_{ref}}}
460 // END_LATEX
461 // where qc is the total charge collected in the current time bin and dx is the length 
462 // of the time bin. For the moment (Jan 20 2009) only pad row cross corrections are 
463 // considered for the charge but none are applied for drift velocity variations along 
464 // the drift region or assymetry of the TRF
465 // 
466 // Author : Alex Bercuci <A.Bercuci@gsi.de>
467 //
468   Float_t dq = 0.;
469   if(fClusters[ic]) dq += TMath::Abs(fClusters[ic]->GetQ());
470   if(fClusters[ic+kNtb]) dq += TMath::Abs(fClusters[ic+kNtb]->GetQ());
471   if(dq<1.e-3 || fdX < 1.e-3) return 0.;
472
473   return dq/fdX/TMath::Sqrt(1. + fYfit[1]*fYfit[1] + fZref[1]*fZref[1]);
474 }
475
476 //____________________________________________________________________
477 Float_t* AliTRDseedV1::GetProbability()
478 {       
479 // Fill probability array for tracklet from the DB.
480 //
481 // Parameters
482 //
483 // Output
484 //   returns pointer to the probability array and 0x0 if missing DB access 
485 //
486 // Detailed description
487
488   
489   // retrive calibration db
490   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
491   if (!calibration) {
492     AliError("No access to calibration data");
493     return 0x0;
494   }
495
496   if (!fReconstructor) {
497     AliError("Reconstructor not set.");
498     return 0x0;
499   }
500
501   // Retrieve the CDB container class with the parametric detector response
502   const AliTRDCalPID *pd = calibration->GetPIDObject(fReconstructor->GetPIDMethod());
503   if (!pd) {
504     AliError("No access to AliTRDCalPID object");
505     return 0x0;
506   }
507   //AliInfo(Form("Method[%d] : %s", fReconstructor->GetRecoParam() ->GetPIDMethod(), pd->IsA()->GetName()));
508
509   // calculate tracklet length TO DO
510   Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
511   /// TMath::Sqrt((1.0 - fSnp[iPlane]*fSnp[iPlane]) / (1.0 + fTgl[iPlane]*fTgl[iPlane]));
512   
513   //calculate dE/dx
514   CookdEdx(fReconstructor->GetNdEdxSlices());
515   
516   // Sets the a priori probabilities
517   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
518     fProb[ispec] = pd->GetProbability(ispec, fMom, &fdEdx[0], length, GetPlane());      
519   }
520
521   return &fProb[0];
522 }
523
524 //____________________________________________________________________
525 Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
526 {
527   //
528   // Returns a quality measurement of the current seed
529   //
530
531   Float_t zcorr = kZcorr ? fTilt * (fZfit[0] - fZref[0]) : 0.;
532   return 
533       .5 * TMath::Abs(18.0 - fN2)
534     + 10.* TMath::Abs(fYfit[1] - fYref[1])
535     + 5. * TMath::Abs(fYfit[0] - fYref[0] + zcorr)
536     + 2. * TMath::Abs(fZfit[0] - fZref[0]) / fPadLength;
537 }
538
539 //____________________________________________________________________
540 void AliTRDseedV1::GetCovAt(Double_t x, Double_t *cov) const
541 {
542 // Computes covariance in the y-z plane at radial point x (in tracking coordinates) 
543 // and returns the results in the preallocated array cov[3] as :
544 //   cov[0] = Var(y)
545 //   cov[1] = Cov(yz)
546 //   cov[2] = Var(z)
547 //
548 // Details
549 //
550 // For the linear transformation
551 // BEGIN_LATEX
552 // Y = T_{x} X^{T}
553 // END_LATEX
554 //   The error propagation has the general form
555 // BEGIN_LATEX
556 // C_{Y} = T_{x} C_{X} T_{x}^{T} 
557 // END_LATEX
558 //  We apply this formula 2 times. First to calculate the covariance of the tracklet 
559 // at point x we consider: 
560 // BEGIN_LATEX
561 // T_{x} = (1 x); X=(y0 dy/dx); C_{X}=#(){#splitline{Var(y0) Cov(y0, dy/dx)}{Cov(y0, dy/dx) Var(dy/dx)}} 
562 // END_LATEX
563 // and secondly to take into account the tilt angle
564 // BEGIN_LATEX
565 // T_{#alpha} = #(){#splitline{cos(#alpha) __ sin(#alpha)}{-sin(#alpha) __ cos(#alpha)}}; X=(y z); C_{X}=#(){#splitline{Var(y)    0}{0   Var(z)}} 
566 // END_LATEX
567 //
568 // using simple trigonometrics one can write for this last case
569 // BEGIN_LATEX
570 // 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})}} 
571 // END_LATEX
572 // which can be aproximated for small alphas (2 deg) with
573 // BEGIN_LATEX
574 // 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}}} 
575 // END_LATEX
576 //
577 // before applying the tilt rotation we also apply systematic uncertainties to the tracklet 
578 // position which can be tunned from outside via the AliTRDrecoParam::SetSysCovMatrix(). They might 
579 // account for extra misalignment/miscalibration uncertainties. 
580 //
581 // Author :
582 // Alex Bercuci <A.Bercuci@gsi.de> 
583 // Date : Jan 8th 2009
584 //
585
586
587   Double_t xr     = fX0-x; 
588   Double_t sy2    = fCov[0] +2.*xr*fCov[1] + xr*xr*fCov[2];
589   Double_t sz2    = fPadLength*fPadLength/12.;
590
591   // insert systematic uncertainties
592   Double_t sys[15];
593   fReconstructor->GetRecoParam()->GetSysCovMatrix(sys);
594   sy2 += sys[0];
595   sz2 += sys[1];
596
597   // rotate covariance matrix
598   Double_t t2 = fTilt*fTilt;
599   Double_t correction = 1./(1. + t2);
600   cov[0] = (sy2+t2*sz2)*correction;
601   cov[1] = fTilt*(sz2 - sy2)*correction;
602   cov[2] = (t2*sy2+sz2)*correction;
603 }
604
605
606 //____________________________________________________________________
607 void AliTRDseedV1::Calibrate()
608 {
609 // Retrieve calibration and position parameters from OCDB. 
610 // The following information are used
611 //  - detector index
612 //  - column and row position of first attached cluster. If no clusters are attached 
613 // to the tracklet a random central chamber position (c=70, r=7) will be used.
614 //
615 // The following information is cached in the tracklet
616 //   t0 (trigger delay)
617 //   drift velocity
618 //   PRF width
619 //   omega*tau = tg(a_L)
620 //   diffusion coefficients (longitudinal and transversal)
621 //
622 // Author :
623 // Alex Bercuci <A.Bercuci@gsi.de> 
624 // Date : Jan 8th 2009
625 //
626
627   AliCDBManager *cdb = AliCDBManager::Instance();
628   if(cdb->GetRun() < 0){
629     AliError("OCDB manager not properly initialized");
630     return;
631   }
632
633   AliTRDcalibDB *calib = AliTRDcalibDB::Instance();
634   AliTRDCalROC  *vdROC = calib->GetVdriftROC(fDet),
635                 *t0ROC = calib->GetT0ROC(fDet);;
636   const AliTRDCalDet *vdDet = calib->GetVdriftDet();
637   const AliTRDCalDet *t0Det = calib->GetT0Det();
638
639   Int_t col = 70, row = 7;
640   AliTRDcluster **c = &fClusters[0];
641   if(fN2){ 
642     Int_t ic = 0;
643     while (ic<kNTimeBins && !(*c)){ic++; c++;} 
644     if(*c){
645       col = (*c)->GetPadCol();
646       row = (*c)->GetPadRow();
647     }
648   }
649
650   fT0    = t0Det->GetValue(fDet) + t0ROC->GetValue(col,row);
651   fVD    = vdDet->GetValue(fDet) * vdROC->GetValue(col, row);
652   fS2PRF = calib->GetPRFWidth(fDet, col, row); fS2PRF *= fS2PRF;
653   fExB   = AliTRDCommonParam::Instance()->GetOmegaTau(fVD);
654   AliTRDCommonParam::Instance()->GetDiffCoeff(fDiffL,
655   fDiffT, fVD);
656   SetBit(kCalib, kTRUE);
657 }
658
659 //____________________________________________________________________
660 void AliTRDseedV1::SetOwner()
661 {
662   //AliInfo(Form("own [%s] fOwner[%s]", own?"YES":"NO", fOwner?"YES":"NO"));
663   
664   if(TestBit(kOwner)) return;
665   for(int ic=0; ic<kNTimeBins; ic++){
666     if(!fClusters[ic]) continue;
667     fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
668   }
669   SetBit(kOwner);
670 }
671
672 //____________________________________________________________________
673 Bool_t  AliTRDseedV1::AttachClustersIter(AliTRDtrackingChamber *chamber, Float_t quality, Bool_t kZcorr, AliTRDcluster *c)
674 {
675   //
676   // Iterative process to register clusters to the seed.
677   // In iteration 0 we try only one pad-row and if quality not
678   // sufficient we try 2 pad-rows (about 5% of tracks cross 2 pad-rows)
679   //
680   // debug level 7
681   //
682   
683   if(!fReconstructor->GetRecoParam() ){
684     AliError("Seed can not be used without a valid RecoParam.");
685     return kFALSE;
686   }
687
688   AliTRDchamberTimeBin *layer = 0x0;
689   if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7){
690     AliTRDtrackingChamber ch(*chamber);
691     ch.SetOwner(); 
692     TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker);
693     cstreamer << "AttachClustersIter"
694       << "chamber.="   << &ch
695       << "tracklet.="  << this
696       << "\n";  
697   }
698
699   Float_t  tquality;
700   Double_t kroady = fReconstructor->GetRecoParam() ->GetRoad1y();
701   Double_t kroadz = fPadLength * .5 + 1.;
702   
703   // initialize configuration parameters
704   Float_t zcorr = kZcorr ? fTilt * (fZfit[0] - fZref[0]) : 0.;
705   Int_t   niter = kZcorr ? 1 : 2;
706   
707   Double_t yexp, zexp;
708   Int_t ncl = 0;
709   // start seed update
710   for (Int_t iter = 0; iter < niter; iter++) {
711     ncl = 0;
712     for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
713       if(!(layer = chamber->GetTB(iTime))) continue;
714       if(!Int_t(*layer)) continue;
715       
716       // define searching configuration
717       Double_t dxlayer = layer->GetX() - fX0;
718       if(c){
719         zexp = c->GetZ();
720         //Try 2 pad-rows in second iteration
721         if (iter > 0) {
722           zexp = fZref[0] + fZref[1] * dxlayer - zcorr;
723           if (zexp > c->GetZ()) zexp = c->GetZ() + fPadLength*0.5;
724           if (zexp < c->GetZ()) zexp = c->GetZ() - fPadLength*0.5;
725         }
726       } else zexp = fZref[0] + (kZcorr ? fZref[1] * dxlayer : 0.);
727       yexp  = fYref[0] + fYref[1] * dxlayer - zcorr;
728       
729       // Get and register cluster
730       Int_t    index = layer->SearchNearestCluster(yexp, zexp, kroady, kroadz);
731       if (index < 0) continue;
732       AliTRDcluster *cl = (*layer)[index];
733       
734       fIndexes[iTime]  = layer->GetGlobalIndex(index);
735       fClusters[iTime] = cl;
736 //       fY[iTime]        = cl->GetY();
737 //       fZ[iTime]        = cl->GetZ();
738       ncl++;
739     }
740     if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d ncl [%d] = %d", iter, fDet, ncl));
741     
742     if(ncl>1){  
743       // calculate length of the time bin (calibration aware)
744       Int_t irp = 0; Float_t x[2]={0., 0.}; Int_t tb[2] = {0,0};
745       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
746         if(!fClusters[iTime]) continue;
747         x[irp]  = fClusters[iTime]->GetX();
748         tb[irp] = iTime;
749         irp++;
750         if(irp==2) break;
751       } 
752       Int_t dtb = tb[1] - tb[0];
753       fdX = dtb ? (x[0] - x[1]) / dtb : 0.15;
754
755       // update X0 from the clusters (calibration/alignment aware)
756       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
757         if(!(layer = chamber->GetTB(iTime))) continue;
758         if(!layer->IsT0()) continue;
759         if(fClusters[iTime]){ 
760           fX0 = fClusters[iTime]->GetX();
761           break;
762         } else { // we have to infere the position of the anode wire from the other clusters
763           for (Int_t jTime = iTime+1; jTime < AliTRDtrackerV1::GetNTimeBins(); jTime++) {
764             if(!fClusters[jTime]) continue;
765             fX0 = fClusters[jTime]->GetX() + fdX * (jTime - iTime);
766             break;
767           }
768         }
769       } 
770       
771       // update YZ reference point
772       // TODO
773       
774       // update x reference positions (calibration/alignment aware)
775 //       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
776 //         if(!fClusters[iTime]) continue;
777 //         fX[iTime] = fX0 - fClusters[iTime]->GetX();
778 //       } 
779       
780       FitMI();
781     }
782     if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d nclFit [%d] = %d", iter, fDet, fN2));
783     
784     if(IsOK()){
785       tquality = GetQuality(kZcorr);
786       if(tquality < quality) break;
787       else quality = tquality;
788     }
789     kroadz *= 2.;
790   } // Loop: iter
791   if (!IsOK()) return kFALSE;
792
793   if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=1) CookLabels();
794
795   // load calibration params
796   Calibrate();
797   UpdateUsed();
798   return kTRUE; 
799 }
800
801 //____________________________________________________________________
802 Bool_t  AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *chamber, Bool_t tilt)
803 {
804   //
805   // Projective algorithm to attach clusters to seeding tracklets
806   //
807   // Parameters
808   //
809   // Output
810   //
811   // Detailed description
812   // 1. Collapse x coordinate for the full detector plane
813   // 2. truncated mean on y (r-phi) direction
814   // 3. purge clusters
815   // 4. truncated mean on z direction
816   // 5. purge clusters
817   // 6. fit tracklet
818   //    
819   Bool_t kPRINT = kFALSE;
820   if(!fReconstructor->GetRecoParam() ){
821     AliError("Seed can not be used without a valid RecoParam.");
822     return kFALSE;
823   }
824   // Initialize reco params for this tracklet
825   // 1. first time bin in the drift region
826   Int_t t0 = 4;
827   Int_t kClmin = Int_t(fReconstructor->GetRecoParam() ->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
828
829   Double_t syRef  = TMath::Sqrt(fRefCov[0]);
830   //define roads
831   Double_t kroady = 1.; 
832   //fReconstructor->GetRecoParam() ->GetRoad1y();
833   Double_t kroadz = fPadLength * 1.5 + 1.;
834   if(kPRINT) printf("AttachClusters() sy[%f] road[%f]\n", syRef, kroady);
835
836   // working variables
837   const Int_t kNrows = 16;
838   AliTRDcluster *clst[kNrows][kNTimeBins];
839   Double_t cond[4], dx, dy, yt, zt,
840     yres[kNrows][kNTimeBins];
841   Int_t idxs[kNrows][kNTimeBins], ncl[kNrows], ncls = 0;
842   memset(ncl, 0, kNrows*sizeof(Int_t));
843   memset(clst, 0, kNrows*kNTimeBins*sizeof(AliTRDcluster*));
844
845   // Do cluster projection
846   AliTRDcluster *c = 0x0;
847   AliTRDchamberTimeBin *layer = 0x0;
848   Bool_t kBUFFER = kFALSE;
849   for (Int_t it = 0; it < AliTRDtrackerV1::GetNTimeBins(); it++) {
850     if(!(layer = chamber->GetTB(it))) continue;
851     if(!Int_t(*layer)) continue;
852     
853     dx   = fX0 - layer->GetX();
854     yt = fYref[0] - fYref[1] * dx;
855     zt = fZref[0] - fZref[1] * dx;
856     if(kPRINT) printf("\t%2d dx[%f] yt[%f] zt[%f]\n", it, dx, yt, zt);
857
858     // select clusters on a 5 sigmaKalman level
859     cond[0] = yt; cond[2] = kroady;
860     cond[1] = zt; cond[3] = kroadz;
861     Int_t n=0, idx[6];
862     layer->GetClusters(cond, idx, n, 6);
863     for(Int_t ic = n; ic--;){
864       c  = (*layer)[idx[ic]];
865       dy = yt - c->GetY();
866       dy += tilt ? fTilt * (c->GetZ() - zt) : 0.;
867       // select clusters on a 3 sigmaKalman level
868 /*      if(tilt && TMath::Abs(dy) > 3.*syRef){ 
869         printf("too large !!!\n");
870         continue;
871       }*/
872       Int_t r = c->GetPadRow();
873       if(kPRINT) printf("\t\t%d dy[%f] yc[%f] r[%d]\n", ic, TMath::Abs(dy), c->GetY(), r);
874       clst[r][ncl[r]] = c;
875       idxs[r][ncl[r]] = idx[ic];
876       yres[r][ncl[r]] = dy;
877       ncl[r]++; ncls++;
878
879       if(ncl[r] >= kNTimeBins) {
880         AliWarning(Form("Cluster candidates reached limit %d. Some may be lost.", kNTimeBins));
881         kBUFFER = kTRUE;
882         break;
883       }
884     }
885     if(kBUFFER) break;
886   }
887   if(kPRINT) printf("Found %d clusters\n", ncls);
888   if(ncls<kClmin) return kFALSE;
889  
890   // analyze each row individualy
891   Double_t mean, syDis;
892   Int_t nrow[] = {0, 0, 0}, nr = 0, lr=-1;
893   for(Int_t ir=kNrows; ir--;){
894     if(!(ncl[ir])) continue;
895     if(lr>0 && lr-ir != 1){
896       if(kPRINT) printf("W - gap in rows attached !!\n"); 
897     }
898     if(kPRINT) printf("\tir[%d] lr[%d] n[%d]\n", ir, lr, ncl[ir]);
899     // Evaluate truncated mean on the y direction
900     if(ncl[ir] > 3) AliMathBase::EvaluateUni(ncl[ir], yres[ir], mean, syDis, Int_t(ncl[ir]*.8));
901     else {
902       mean = 0.; syDis = 0.;
903     } 
904
905     // TODO check mean and sigma agains cluster resolution !!
906     if(kPRINT) printf("\tr[%2d] m[%f %5.3fsigma] s[%f]\n", ir, mean, TMath::Abs(mean/syRef), syDis);
907     // select clusters on a 3 sigmaDistr level
908     Bool_t kFOUND = kFALSE;
909     for(Int_t ic = ncl[ir]; ic--;){
910       if(yres[ir][ic] - mean > 3. * syDis){ 
911         clst[ir][ic] = 0x0; continue;
912       }
913       nrow[nr]++; kFOUND = kTRUE;
914     }
915     // exit loop
916     if(kFOUND) nr++; 
917     lr = ir; if(nr>=3) break;
918   }
919   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]);
920
921   // classify cluster rows
922   Int_t row = -1;
923   switch(nr){
924   case 1:
925     row = lr;
926     break;
927   case 2:
928     SetBit(kRowCross, kTRUE); // mark pad row crossing
929     if(nrow[0] > nrow[1]){ row = lr+1; lr = -1;}
930     else{ 
931       row = lr; lr = 1;
932       nrow[2] = nrow[1];
933       nrow[1] = nrow[0];
934       nrow[0] = nrow[2];
935     }
936     break;
937   case 3:
938     SetBit(kRowCross, kTRUE); // mark pad row crossing
939     break;
940   }
941   if(kPRINT) printf("\trow[%d] n[%d]\n\n", row, nrow[0]);
942   if(row<0) return kFALSE;
943
944   // Select and store clusters 
945   // We should consider here :
946   //  1. How far is the chamber boundary
947   //  2. How big is the mean
948   fN2 = 0;
949   for (Int_t ir = 0; ir < nr; ir++) {
950     Int_t jr = row + ir*lr; 
951     if(kPRINT) printf("\tattach %d clusters for row %d\n", ncl[jr], jr);
952     for (Int_t ic = 0; ic < ncl[jr]; ic++) {
953       if(!(c = clst[jr][ic])) continue;
954       Int_t it = c->GetPadTime();
955       // TODO proper indexing of clusters !!
956       fIndexes[it+kNtb*ir]  = chamber->GetTB(it)->GetGlobalIndex(idxs[jr][ic]);
957       fClusters[it+kNtb*ir] = c;
958   
959       //printf("\tid[%2d] it[%d] idx[%d]\n", ic, it, fIndexes[it]);
960   
961       fN2++;
962     }
963   }  
964
965   // number of minimum numbers of clusters expected for the tracklet
966   if (fN2 < kClmin){
967     AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin));
968     fN2 = 0;
969     return kFALSE;
970   }
971
972   // update used clusters and select
973   fNUsed = 0;
974   for (Int_t it = 0; it < AliTRDtrackerV1::GetNTimeBins(); it++) {
975     if(fClusters[it] && fClusters[it]->IsUsed()) fNUsed++;
976     if(fClusters[it+kNtb] && fClusters[it+kNtb]->IsUsed()) fNUsed++;
977   }
978   if (fN2-fNUsed < kClmin){
979     //AliWarning(Form("Too many clusters already in use %d (from %d).", fNUsed, fN2));
980     fN2 = 0;
981     return kFALSE;
982   }
983
984   // Load calibration parameters for this tracklet  
985   Calibrate();
986
987   // calculate dx for time bins in the drift region (calibration aware)
988   Int_t irp = 0; Float_t x[2] = {0.,0.}; Int_t tb[2]={0,0};
989   for (Int_t it = t0; it < AliTRDtrackerV1::GetNTimeBins(); it++) {
990     if(!fClusters[it]) continue;
991     x[irp]  = fClusters[it]->GetX();
992     tb[irp] = it;
993     irp++;
994     if(irp==2) break;
995   }  
996   Int_t dtb = tb[1] - tb[0];
997   fdX = dtb ? (x[0] - x[1]) / dtb : 0.15;
998
999   // update X0 from the clusters (calibration/alignment aware) TODO remove dependence on x0 !!
1000   for (Int_t it = 0; it < AliTRDtrackerV1::GetNTimeBins(); it++) {
1001     if(!(layer = chamber->GetTB(it))) continue;
1002     if(!layer->IsT0()) continue;
1003     if(fClusters[it]){ 
1004       fX0 = fClusters[it]->GetX();
1005       break;
1006     } else { // we have to infere the position of the anode wire from the other clusters
1007       for (Int_t jt = it+1; jt < AliTRDtrackerV1::GetNTimeBins(); jt++) {
1008         if(!fClusters[jt]) continue;
1009         fX0 = fClusters[jt]->GetX() + fdX * (jt - it);
1010         break;
1011       }
1012     }
1013   }     
1014
1015   return kTRUE;
1016 }
1017
1018 //____________________________________________________________
1019 void AliTRDseedV1::Bootstrap(const AliTRDReconstructor *rec)
1020 {
1021 //   Fill in all derived information. It has to be called after recovery from file or HLT.
1022 //   The primitive data are
1023 //   - list of clusters
1024 //   - detector (as the detector will be removed from clusters)
1025 //   - position of anode wire (fX0) - temporary
1026 //   - track reference position and direction
1027 //   - momentum of the track
1028 //   - time bin length [cm]
1029 // 
1030 //   A.Bercuci <A.Bercuci@gsi.de> Oct 30th 2008
1031 //
1032   fReconstructor = rec;
1033   AliTRDgeometry g;
1034   AliTRDpadPlane *pp = g.GetPadPlane(fDet);
1035   fTilt      = TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle());
1036   fPadLength = pp->GetLengthIPad();
1037   //fSnp = fYref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);
1038   //fTgl = fZref[1];
1039   fN2 = 0;// fMPads = 0.;
1040   AliTRDcluster **cit = &fClusters[0];
1041   for(Int_t ic = kNTimeBins; ic--; cit++){
1042     if(!(*cit)) return;
1043     fN2++;
1044 /*    fX[ic] = (*cit)->GetX() - fX0;
1045     fY[ic] = (*cit)->GetY();
1046     fZ[ic] = (*cit)->GetZ();*/
1047   }
1048   //Update(); // 
1049   Fit();
1050   CookLabels();
1051   GetProbability();
1052 }
1053
1054
1055 //____________________________________________________________________
1056 Bool_t AliTRDseedV1::Fit(Bool_t tilt, Int_t errors)
1057 {
1058   //
1059   // Linear fit of the tracklet
1060   //
1061   // Parameters :
1062   //
1063   // Output :
1064   //  True if successful
1065   //
1066   // Detailed description
1067   // 2. Check if tracklet crosses pad row boundary
1068   // 1. Calculate residuals in the y (r-phi) direction
1069   // 3. Do a Least Square Fit to the data
1070   //
1071
1072   if(!IsCalibrated()){
1073     AliWarning("Tracklet fit failed. Call Calibrate().");
1074     return kFALSE;
1075   }
1076
1077   const Int_t kClmin = 8;
1078
1079
1080   // cluster error parametrization parameters 
1081   // 1. sy total charge
1082   const Float_t sq0inv = 0.019962; // [1/q0]
1083   const Float_t sqb    = 1.0281564;    //[cm]
1084   // 2. sy for the PRF
1085   const Float_t scy[AliTRDgeometry::kNlayer][4] = {
1086     {2.827e-02, 9.600e-04, 4.296e-01, 2.271e-02},
1087     {2.952e-02,-2.198e-04, 4.146e-01, 2.339e-02},
1088     {3.090e-02, 1.514e-03, 4.020e-01, 2.402e-02},
1089     {3.260e-02,-2.037e-03, 3.946e-01, 2.509e-02},
1090     {3.439e-02,-3.601e-04, 3.883e-01, 2.623e-02},
1091     {3.510e-02, 2.066e-03, 3.651e-01, 2.588e-02},
1092   };
1093   // 3. sy parallel to the track
1094   const Float_t sy0 =  2.649e-02; // [cm]
1095   const Float_t sya = -8.864e-04; // [cm]
1096   const Float_t syb = -2.435e-01; // [cm]
1097
1098   // 4. sx parallel to the track
1099   const Float_t sxgc = 5.427e-02;
1100   const Float_t sxgm = 7.783e-01;
1101   const Float_t sxgs = 2.743e-01;
1102   const Float_t sxe0 =-2.065e+00;
1103   const Float_t sxe1 =-2.978e-02;
1104
1105   // 5. sx perpendicular to the track
1106 //   const Float_t sxd0 = 1.881e-02;
1107 //   const Float_t sxd1 =-4.101e-01;
1108 //   const Float_t sxd2 = 1.572e+00;
1109
1110   // get track direction
1111   Double_t y0   = fYref[0];
1112   Double_t dydx = fYref[1]; 
1113   Double_t z0   = fZref[0];
1114   Double_t dzdx = fZref[1];
1115   Double_t yt, zt;
1116
1117   const Int_t kNtb = AliTRDtrackerV1::GetNTimeBins();
1118   // calculation of tg^2(phi - a_L) and tg^2(a_L)
1119   Double_t tgg = (dydx-fExB)/(1.+dydx*fExB); tgg *= tgg;
1120   //Double_t exb2= fExB*fExB;
1121
1122   //AliTRDtrackerV1::AliTRDLeastSquare fitterZ;
1123   TLinearFitter  fitterY(1, "pol1");
1124   // convertion factor from square to gauss distribution for sigma
1125   //Double_t convert = 1./TMath::Sqrt(12.);
1126   
1127   // book cluster information
1128   Double_t qc[kNTimeBins], xc[kNTimeBins], yc[kNTimeBins], zc[kNTimeBins], sy[kNTimeBins];
1129
1130   Int_t ily = AliTRDgeometry::GetLayer(fDet);
1131   Int_t fN = 0;
1132   AliTRDcluster *c=0x0, **jc = &fClusters[0];
1133   for (Int_t ic=0; ic<kNtb; ic++, ++jc) {
1134     //zRow[ic] = -1;
1135     xc[ic]  = -1.;
1136     yc[ic]  = 999.;
1137     zc[ic]  = 999.;
1138     sy[ic]  = 0.;
1139     //sz[ic]  = 0.;
1140     if(!(c = (*jc))) continue;
1141     if(!c->IsInChamber()) continue;
1142
1143     Float_t w = 1.;
1144     if(c->GetNPads()>4) w = .5;
1145     if(c->GetNPads()>5) w = .2;
1146
1147     //zRow[fN] = c->GetPadRow();
1148     qc[fN]   = TMath::Abs(c->GetQ());
1149     // correct cluster position for PRF and v drift
1150     //Int_t jc = TMath::Max(fN-3, 0);
1151     //xc[fN]   = c->GetXloc(fT0, fVD, &qc[jc], &xc[jc]/*, z0 - c->GetX()*dzdx*/);
1152     //Double_t s2 = fS2PRF + fDiffL*fDiffL*xc[fN]/(1.+2.*exb2)+tgg*xc[fN]*xc[fN]*exb2/12.;
1153     //yc[fN]   = c->GetYloc(s2, fPadLength, xc[fN], fExB);
1154     
1155     // uncalibrated cluster correction 
1156     // TODO remove
1157     Double_t x, y; GetClusterXY(c, x, y);
1158     xc[fN]   = fX0 - x;
1159     yc[fN]   = y;
1160     zc[fN]   = c->GetZ();
1161
1162     // extrapolated y value for the track
1163     yt = y0 - xc[fN]*dydx; 
1164     // extrapolated z value for the track
1165     zt = z0 - xc[fN]*dzdx; 
1166     // tilt correction
1167     if(tilt) yc[fN] -= fTilt*(zc[fN] - zt); 
1168
1169     // ELABORATE CLUSTER ERROR
1170     // TODO to be moved to AliTRDcluster
1171     // basic y error (|| to track).
1172     sy[fN]  = xc[fN] < AliTRDgeometry::CamHght() ? 2. : sy0 + sya*TMath::Exp(1./(xc[fN]+syb));
1173     //printf("cluster[%d]\n\tsy[0] = %5.3e [um]\n", fN,  sy[fN]*1.e4);
1174     // y error due to total charge
1175     sy[fN] += sqb*(1./qc[fN] - sq0inv);
1176     //printf("\tsy[1] = %5.3e [um]\n", sy[fN]*1.e4);
1177     // y error due to PRF
1178     sy[fN] += scy[ily][0]*TMath::Gaus(c->GetCenter(), scy[ily][1], scy[ily][2]) - scy[ily][3];
1179     //printf("\tsy[2] = %5.3e [um]\n", sy[fN]*1.e4);
1180
1181     sy[fN] *= sy[fN];
1182
1183     // ADD ERROR ON x
1184     // error of drift length parallel to the track
1185     Double_t sx = sxgc*TMath::Gaus(xc[fN], sxgm, sxgs) + TMath::Exp(sxe0+sxe1*xc[fN]); // [cm]
1186     //printf("\tsx[0] = %5.3e [um]\n", sx*1.e4);
1187     // error of drift length perpendicular to the track
1188     //sx += sxd0 + sxd1*d + sxd2*d*d;
1189     sx *= sx; // square sx
1190
1191     // add error from ExB 
1192     if(errors>0) sy[fN] += fExB*fExB*sx;
1193     //printf("\tsy[3] = %5.3e [um^2]\n", sy[fN]*1.e8);
1194
1195     // global radial error due to misalignment/miscalibration
1196     Double_t sx0  = 0.; sx0 *= sx0;
1197     // add sx contribution to sy due to track angle
1198     if(errors>1) sy[fN] += tgg*(sx+sx0);
1199     // TODO we should add tilt pad correction here
1200     //printf("\tsy[4] = %5.3e [um^2]\n", sy[fN]*1.e8);
1201     c->SetSigmaY2(sy[fN]);
1202
1203     sy[fN]  = TMath::Sqrt(sy[fN]);
1204     fitterY.AddPoint(&xc[fN], yc[fN], sy[fN]);
1205     fN++;
1206   }
1207   // to few clusters
1208   if (fN < kClmin) return kFALSE; 
1209
1210   // fit XY
1211   fitterY.Eval();
1212   fYfit[0] = fitterY.GetParameter(0);
1213   fYfit[1] = -fitterY.GetParameter(1);
1214   // store covariance
1215   Double_t *p = fitterY.GetCovarianceMatrix();
1216   fCov[0] = p[0]; // variance of y0
1217   fCov[1] = p[1]; // covariance of y0, dydx
1218   fCov[2] = p[3]; // variance of dydx
1219   // the ref radial position is set at the minimum of 
1220   // the y variance of the tracklet
1221   fX = -fCov[1]/fCov[2]; //fXref = fX0 - fXref;
1222
1223   // fit XZ
1224   if(IsRowCross()){ 
1225     // TODO pad row cross position estimation !!!
1226     //AliInfo(Form("Padrow cross in detector %d", fDet));
1227     fZfit[0] = .5*(zc[0]+zc[fN-1]); fZfit[1] = 0.;
1228   } else {
1229     fZfit[0] = zc[0]; fZfit[1] = 0.;
1230   }
1231
1232
1233 //   // determine z offset of the fit
1234 //   Float_t zslope = 0.;
1235 //   Int_t nchanges = 0, nCross = 0;
1236 //   if(nz==2){ // tracklet is crossing pad row
1237 //     // Find the break time allowing one chage on pad-rows
1238 //     // with maximal number of accepted clusters
1239 //     Int_t padRef = zRow[0];
1240 //     for (Int_t ic=1; ic<fN; ic++) {
1241 //       if(zRow[ic] == padRef) continue;
1242 //       
1243 //       // debug
1244 //       if(zRow[ic-1] == zRow[ic]){
1245 //         printf("ERROR in pad row change!!!\n");
1246 //       }
1247 //     
1248 //       // evaluate parameters of the crossing point
1249 //       Float_t sx = (xc[ic-1] - xc[ic])*convert;
1250 //       fCross[0] = .5 * (xc[ic-1] + xc[ic]);
1251 //       fCross[2] = .5 * (zc[ic-1] + zc[ic]);
1252 //       fCross[3] = TMath::Max(dzdx * sx, .01);
1253 //       zslope    = zc[ic-1] > zc[ic] ? 1. : -1.;
1254 //       padRef    = zRow[ic];
1255 //       nCross    = ic;
1256 //       nchanges++;
1257 //     }
1258 //   }
1259 // 
1260 //   // condition on nCross and reset nchanges TODO
1261 // 
1262 //   if(nchanges==1){
1263 //     if(dzdx * zslope < 0.){
1264 //       AliInfo("Tracklet-Track mismatch in dzdx. TODO.");
1265 //     }
1266 // 
1267 // 
1268 //     //zc[nc] = fitterZ.GetFunctionParameter(0); 
1269 //     fCross[1] = fYfit[0] - fCross[0] * fYfit[1];
1270 //     fCross[0] = fX0 - fCross[0];
1271 //   }
1272
1273   UpdateUsed();
1274   return kTRUE;
1275 }
1276
1277
1278
1279
1280
1281
1282 //_____________________________________________________________________________
1283 void AliTRDseedV1::FitMI()
1284 {
1285 //
1286 // Fit the seed.
1287 // Marian Ivanov's version 
1288 //
1289 // linear fit on the y direction with respect to the reference direction. 
1290 // The residuals for each x (x = xc - x0) are deduced from:
1291 // dy = y - yt             (1)
1292 // the tilting correction is written :
1293 // y = yc + h*(zc-zt)      (2)
1294 // yt = y0+dy/dx*x         (3)
1295 // zt = z0+dz/dx*x         (4)
1296 // from (1),(2),(3) and (4)
1297 // dy = yc - y0 - (dy/dx + h*dz/dx)*x + h*(zc-z0)
1298 // the last term introduces the correction on y direction due to tilting pads. There are 2 ways to account for this:
1299 // 1. use tilting correction for calculating the y
1300 // 2. neglect tilting correction here and account for it in the error parametrization of the tracklet.
1301   const Float_t kRatio  = 0.8;
1302   const Int_t   kClmin  = 5;
1303   const Float_t kmaxtan = 2;
1304
1305   if (TMath::Abs(fYref[1]) > kmaxtan){
1306                 //printf("Exit: Abs(fYref[1]) = %3.3f, kmaxtan = %3.3f\n", TMath::Abs(fYref[1]), kmaxtan);
1307                 return;              // Track inclined too much
1308         }
1309
1310   Float_t  sigmaexp  = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction
1311   Float_t  ycrosscor = fPadLength * fTilt * 0.5;           // Y correction for crossing 
1312   Int_t fNChange = 0;
1313
1314   Double_t sumw;
1315   Double_t sumwx;
1316   Double_t sumwx2;
1317   Double_t sumwy;
1318   Double_t sumwxy;
1319   Double_t sumwz;
1320   Double_t sumwxz;
1321
1322         // Buffering: Leave it constant fot Performance issues
1323   Int_t    zints[kNtb];            // Histograming of the z coordinate 
1324                                          // Get 1 and second max probable coodinates in z
1325   Int_t    zouts[2*kNtb];       
1326   Float_t  allowedz[kNtb];         // Allowed z for given time bin
1327   Float_t  yres[kNtb];             // Residuals from reference
1328   //Float_t  anglecor = fTilt * fZref[1];  // Correction to the angle
1329   
1330   Float_t pos[3*kNtb]; memset(pos, 0, 3*kNtb*sizeof(Float_t));
1331   Float_t *fX = &pos[0], *fY = &pos[kNtb], *fZ = &pos[2*kNtb];
1332   
1333   Int_t fN  = 0; AliTRDcluster *c = 0x0; 
1334   fN2 = 0;
1335   for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1336     yres[i] = 10000.0;
1337     if (!(c = fClusters[i])) continue;
1338     if(!c->IsInChamber()) continue;
1339     // Residual y
1340     //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] + fTilt*(fZ[i] - fZref[0]);
1341     fX[i] = fX0 - c->GetX();
1342     fY[i] = c->GetY();
1343     fZ[i] = c->GetZ();
1344     yres[i] = fY[i] - fTilt*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
1345     zints[fN] = Int_t(fZ[i]);
1346     fN++;
1347   }
1348
1349   if (fN < kClmin){
1350     //printf("Exit fN < kClmin: fN = %d\n", fN);
1351     return; 
1352   }
1353   Int_t nz = AliTRDtrackerV1::Freq(fN, zints, zouts, kFALSE);
1354   Float_t fZProb   = zouts[0];
1355   if (nz <= 1) zouts[3] = 0;
1356   if (zouts[1] + zouts[3] < kClmin) {
1357     //printf("Exit zouts[1] = %d, zouts[3] = %d\n",zouts[1],zouts[3]);
1358     return;
1359   }
1360   
1361   // Z distance bigger than pad - length
1362   if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) zouts[3] = 0;
1363   
1364   Int_t  breaktime = -1;
1365   Bool_t mbefore   = kFALSE;
1366   Int_t  cumul[kNtb][2];
1367   Int_t  counts[2] = { 0, 0 };
1368   
1369   if (zouts[3] >= 3) {
1370
1371     //
1372     // Find the break time allowing one chage on pad-rows
1373     // with maximal number of accepted clusters
1374     //
1375     fNChange = 1;
1376     for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1377       cumul[i][0] = counts[0];
1378       cumul[i][1] = counts[1];
1379       if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++;
1380       if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++;
1381     }
1382     Int_t  maxcount = 0;
1383     for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1384       Int_t after  = cumul[AliTRDtrackerV1::GetNTimeBins()][0] - cumul[i][0];
1385       Int_t before = cumul[i][1];
1386       if (after + before > maxcount) { 
1387         maxcount  = after + before; 
1388         breaktime = i;
1389         mbefore   = kFALSE;
1390       }
1391       after  = cumul[AliTRDtrackerV1::GetNTimeBins()-1][1] - cumul[i][1];
1392       before = cumul[i][0];
1393       if (after + before > maxcount) { 
1394         maxcount  = after + before; 
1395         breaktime = i;
1396         mbefore   = kTRUE;
1397       }
1398     }
1399     breaktime -= 1;
1400   }
1401
1402   for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1403     if (i >  breaktime) allowedz[i] =   mbefore  ? zouts[2] : zouts[0];
1404     if (i <= breaktime) allowedz[i] = (!mbefore) ? zouts[2] : zouts[0];
1405   }  
1406
1407   if (((allowedz[0] > allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] < 0)) ||
1408       ((allowedz[0] < allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] > 0))) {
1409     //
1410     // Tracklet z-direction not in correspondance with track z direction 
1411     //
1412     fNChange = 0;
1413     for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1414       allowedz[i] = zouts[0];  // Only longest taken
1415     } 
1416   }
1417   
1418   if (fNChange > 0) {
1419     //
1420     // Cross pad -row tracklet  - take the step change into account
1421     //
1422     for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1423       if (!fClusters[i]) continue; 
1424       if(!fClusters[i]->IsInChamber()) continue;
1425       if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
1426       // Residual y
1427       //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] /*+ fTilt*(fZ[i] - fZref[0])*/;   
1428       yres[i] = fY[i] - fTilt*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
1429 /*      if (TMath::Abs(fZ[i] - fZProb) > 2) {
1430         if (fZ[i] > fZProb) yres[i] += fTilt * fPadLength;
1431         if (fZ[i] < fZProb) yres[i] -= fTilt * fPadLength;
1432       }*/
1433     }
1434   }
1435   
1436   Double_t yres2[kNtb];
1437   Double_t mean;
1438   Double_t sigma;
1439   for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1440     if (!fClusters[i]) continue;
1441     if(!fClusters[i]->IsInChamber()) continue;
1442     if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
1443     yres2[fN2] = yres[i];
1444     fN2++;
1445   }
1446   if (fN2 < kClmin) {
1447                 //printf("Exit fN2 < kClmin: fN2 = %d\n", fN2);
1448     fN2 = 0;
1449     return;
1450   }
1451   AliMathBase::EvaluateUni(fN2,yres2,mean,sigma, Int_t(fN2*kRatio-2.));
1452   if (sigma < sigmaexp * 0.8) {
1453     sigma = sigmaexp;
1454   }
1455   //Float_t fSigmaY = sigma;
1456
1457   // Reset sums
1458   sumw   = 0; 
1459   sumwx  = 0; 
1460   sumwx2 = 0;
1461   sumwy  = 0; 
1462   sumwxy = 0; 
1463   sumwz  = 0;
1464   sumwxz = 0;
1465
1466   fN2    = 0;
1467   Float_t fMeanz = 0;
1468   Float_t fMPads = 0;
1469   fUsable = 0;
1470   for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1471     if (!fClusters[i]) continue;
1472     if (!fClusters[i]->IsInChamber()) continue;
1473     if (TMath::Abs(fZ[i] - allowedz[i]) > 2){fClusters[i] = 0x0; continue;}
1474     if (TMath::Abs(yres[i] - mean) > 4.0 * sigma){fClusters[i] = 0x0;  continue;}
1475     SETBIT(fUsable,i);
1476     fN2++;
1477     fMPads += fClusters[i]->GetNPads();
1478     Float_t weight = 1.0;
1479     if (fClusters[i]->GetNPads() > 4) weight = 0.5;
1480     if (fClusters[i]->GetNPads() > 5) weight = 0.2;
1481    
1482         
1483     Double_t x = fX[i];
1484     //printf("x = %7.3f dy = %7.3f fit %7.3f\n", x, yres[i], fY[i]-yres[i]);
1485     
1486     sumw   += weight; 
1487     sumwx  += x * weight; 
1488     sumwx2 += x*x * weight;
1489     sumwy  += weight * yres[i];  
1490     sumwxy += weight * (yres[i]) * x;
1491     sumwz  += weight * fZ[i];    
1492     sumwxz += weight * fZ[i] * x;
1493
1494   }
1495
1496   if (fN2 < kClmin){
1497                 //printf("Exit fN2 < kClmin(2): fN2 = %d\n",fN2);
1498     fN2 = 0;
1499     return;
1500   }
1501   fMeanz = sumwz / sumw;
1502   Float_t correction = 0;
1503   if (fNChange > 0) {
1504     // Tracklet on boundary
1505     if (fMeanz < fZProb) correction =  ycrosscor;
1506     if (fMeanz > fZProb) correction = -ycrosscor;
1507   }
1508
1509   Double_t det = sumw * sumwx2 - sumwx * sumwx;
1510   fYfit[0]    = (sumwx2 * sumwy  - sumwx * sumwxy) / det;
1511   fYfit[1]    = (sumw   * sumwxy - sumwx * sumwy)  / det;
1512   
1513   fS2Y = 0;
1514   for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1515     if (!TESTBIT(fUsable,i)) continue;
1516     Float_t delta = yres[i] - fYfit[0] - fYfit[1] * fX[i];
1517     fS2Y += delta*delta;
1518   }
1519   fS2Y = TMath::Sqrt(fS2Y / Float_t(fN2-2));
1520         // TEMPORARY UNTIL covariance properly calculated
1521         fS2Y = TMath::Max(fS2Y, Float_t(.1));
1522   
1523   fZfit[0]   = (sumwx2 * sumwz  - sumwx * sumwxz) / det;
1524   fZfit[1]   = (sumw   * sumwxz - sumwx * sumwz)  / det;
1525 //   fYfitR[0] += fYref[0] + correction;
1526 //   fYfitR[1] += fYref[1];
1527 //  fYfit[0]   = fYfitR[0];
1528   fYfit[1]   = -fYfit[1];
1529
1530   UpdateUsed();
1531 }
1532
1533
1534 //___________________________________________________________________
1535 void AliTRDseedV1::Print(Option_t *o) const
1536 {
1537   //
1538   // Printing the seedstatus
1539   //
1540
1541   AliInfo(Form("Det[%3d] Tilt[%+6.2f] Pad[%5.2f]", fDet, fTilt, fPadLength));
1542   AliInfo(Form("N[%2d] Nuse[%2d]", fN2, fNUsed));
1543   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]));
1544   AliInfo(Form("Ref        y[%7.2f] z[%7.2f] dydx[%5.2f] dzdx[%5.2f]", fYref[0], fZref[0], fYref[1], fZref[1]))
1545
1546
1547   if(strcmp(o, "a")!=0) return;
1548
1549   AliTRDcluster* const* jc = &fClusters[0];
1550   for(int ic=0; ic<kNTimeBins; ic++, jc++) {
1551     if(!(*jc)) continue;
1552     (*jc)->Print(o);
1553   }
1554 }
1555
1556
1557 //___________________________________________________________________
1558 Bool_t AliTRDseedV1::IsEqual(const TObject *o) const
1559 {
1560   // Checks if current instance of the class has the same essential members
1561   // as the given one
1562
1563   if(!o) return kFALSE;
1564   const AliTRDseedV1 *inTracklet = dynamic_cast<const AliTRDseedV1*>(o);
1565   if(!inTracklet) return kFALSE;
1566
1567   for (Int_t i = 0; i < 2; i++){
1568     if ( fYref[i] != inTracklet->fYref[i] ) return kFALSE;
1569     if ( fZref[i] != inTracklet->fZref[i] ) return kFALSE;
1570   }
1571   
1572   if ( fS2Y != inTracklet->fS2Y ) return kFALSE;
1573   if ( fTilt != inTracklet->fTilt ) return kFALSE;
1574   if ( fPadLength != inTracklet->fPadLength ) return kFALSE;
1575   
1576   for (Int_t i = 0; i < kNTimeBins; i++){
1577 //     if ( fX[i] != inTracklet->GetX(i) ) return kFALSE;
1578 //     if ( fY[i] != inTracklet->GetY(i) ) return kFALSE;
1579 //     if ( fZ[i] != inTracklet->GetZ(i) ) return kFALSE;
1580     if ( fIndexes[i] != inTracklet->fIndexes[i] ) return kFALSE;
1581   }
1582   if ( fUsable != inTracklet->fUsable ) return kFALSE;
1583
1584   for (Int_t i=0; i < 2; i++){
1585     if ( fYfit[i] != inTracklet->fYfit[i] ) return kFALSE;
1586     if ( fZfit[i] != inTracklet->fZfit[i] ) return kFALSE;
1587     if ( fLabels[i] != inTracklet->fLabels[i] ) return kFALSE;
1588   }
1589   
1590 /*  if ( fMeanz != inTracklet->GetMeanz() ) return kFALSE;
1591   if ( fZProb != inTracklet->GetZProb() ) return kFALSE;*/
1592   if ( fN2 != inTracklet->fN2 ) return kFALSE;
1593   if ( fNUsed != inTracklet->fNUsed ) return kFALSE;
1594   //if ( fFreq != inTracklet->GetFreq() ) return kFALSE;
1595   //if ( fNChange != inTracklet->GetNChange() ) return kFALSE;
1596    
1597   if ( fC != inTracklet->fC ) return kFALSE;
1598   //if ( fCC != inTracklet->GetCC() ) return kFALSE;
1599   if ( fChi2 != inTracklet->fChi2 ) return kFALSE;
1600   //  if ( fChi2Z != inTracklet->GetChi2Z() ) return kFALSE;
1601
1602   if ( fDet != inTracklet->fDet ) return kFALSE;
1603   if ( fMom != inTracklet->fMom ) return kFALSE;
1604   if ( fdX != inTracklet->fdX ) return kFALSE;
1605   
1606   for (Int_t iCluster = 0; iCluster < kNTimeBins; iCluster++){
1607     AliTRDcluster *curCluster = fClusters[iCluster];
1608     AliTRDcluster *inCluster = inTracklet->fClusters[iCluster];
1609     if (curCluster && inCluster){
1610       if (! curCluster->IsEqual(inCluster) ) {
1611         curCluster->Print();
1612         inCluster->Print();
1613         return kFALSE;
1614       }
1615     } else {
1616       // if one cluster exists, and corresponding 
1617       // in other tracklet doesn't - return kFALSE
1618       if(curCluster || inCluster) return kFALSE;
1619     }
1620   }
1621   return kTRUE;
1622 }