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