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