Coverity fixes by Ivana (again, thanks so much!!) implemented
[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"
eb38ed55 39#include <TTreeStream.h>
e4f2f73d 40
41#include "AliLog.h"
42#include "AliMathBase.h"
4ecadb52 43#include "AliRieman.h"
d937ad7a 44#include "AliCDBManager.h"
e4f2f73d 45
4ecadb52 46#include "AliTRDReconstructor.h"
03cef9b2 47#include "AliTRDpadPlane.h"
e4f2f73d 48#include "AliTRDcluster.h"
f3d3af1b 49#include "AliTRDseedV1.h"
50#include "AliTRDtrackV1.h"
e4f2f73d 51#include "AliTRDcalibDB.h"
eb38ed55 52#include "AliTRDchamberTimeBin.h"
53#include "AliTRDtrackingChamber.h"
54#include "AliTRDtrackerV1.h"
e4f2f73d 55#include "AliTRDrecoParam.h"
a076fc2f 56#include "AliTRDCommonParam.h"
d937ad7a 57
0906e73e 58#include "Cal/AliTRDCalPID.h"
d937ad7a 59#include "Cal/AliTRDCalROC.h"
60#include "Cal/AliTRDCalDet.h"
e4f2f73d 61
4ecadb52 62class AliTracker;
63
e4f2f73d 64ClassImp(AliTRDseedV1)
65
66//____________________________________________________________________
ae4e8b84 67AliTRDseedV1::AliTRDseedV1(Int_t det)
3e778975 68 :AliTRDtrackletBase()
4d6aee34 69 ,fkReconstructor(NULL)
70 ,fClusterIter(NULL)
e3cf3d02 71 ,fExB(0.)
72 ,fVD(0.)
73 ,fT0(0.)
74 ,fS2PRF(0.)
75 ,fDiffL(0.)
76 ,fDiffT(0.)
ae4e8b84 77 ,fClusterIdx(0)
7c3eecb8 78 ,fErrorMsg(0)
3e778975 79 ,fN(0)
ae4e8b84 80 ,fDet(det)
b25a5e9e 81 ,fPt(0.)
bcb6fb78 82 ,fdX(0.)
e3cf3d02 83 ,fX0(0.)
84 ,fX(0.)
85 ,fY(0.)
86 ,fZ(0.)
87 ,fS2Y(0.)
88 ,fS2Z(0.)
e3cf3d02 89 ,fChi2(0.)
e4f2f73d 90{
91 //
92 // Constructor
93 //
f301a656 94 memset(fIndexes,0xFF,kNclusters*sizeof(fIndexes[0]));
8d2bec9e 95 memset(fClusters, 0, kNclusters*sizeof(AliTRDcluster*));
2eb10c34 96 memset(fPad, 0, 4*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
16cca13f 105 memset(fRefCov, 0, 7*sizeof(Double_t));
68f9b6bd 106 // stand alone curvature
107 fC[0] = 0.; fC[1] = 0.;
d937ad7a 108 // covariance matrix [diagonal]
109 // default sy = 200um and sz = 2.3 cm
110 fCov[0] = 4.e-4; fCov[1] = 0.; fCov[2] = 5.3;
f29f13a6 111 SetStandAlone(kFALSE);
e4f2f73d 112}
113
114//____________________________________________________________________
0906e73e 115AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
3e778975 116 :AliTRDtrackletBase((AliTRDtrackletBase&)ref)
4d6aee34 117 ,fkReconstructor(NULL)
118 ,fClusterIter(NULL)
e3cf3d02 119 ,fExB(0.)
120 ,fVD(0.)
121 ,fT0(0.)
122 ,fS2PRF(0.)
123 ,fDiffL(0.)
124 ,fDiffT(0.)
ae4e8b84 125 ,fClusterIdx(0)
7c3eecb8 126 ,fErrorMsg(0)
3e778975 127 ,fN(0)
e3cf3d02 128 ,fDet(-1)
b25a5e9e 129 ,fPt(0.)
e3cf3d02 130 ,fdX(0.)
131 ,fX0(0.)
132 ,fX(0.)
133 ,fY(0.)
134 ,fZ(0.)
135 ,fS2Y(0.)
136 ,fS2Z(0.)
e3cf3d02 137 ,fChi2(0.)
e4f2f73d 138{
139 //
140 // Copy Constructor performing a deep copy
141 //
e3cf3d02 142 if(this != &ref){
143 ref.Copy(*this);
144 }
29b87567 145 SetBit(kOwner, kFALSE);
f29f13a6 146 SetStandAlone(ref.IsStandAlone());
fbb2ea06 147}
d9950a5a 148
0906e73e 149
e4f2f73d 150//____________________________________________________________________
151AliTRDseedV1& AliTRDseedV1::operator=(const AliTRDseedV1 &ref)
152{
153 //
154 // Assignment Operator using the copy function
155 //
156
29b87567 157 if(this != &ref){
158 ref.Copy(*this);
159 }
221ab7e0 160 SetBit(kOwner, kFALSE);
161
29b87567 162 return *this;
e4f2f73d 163}
164
165//____________________________________________________________________
166AliTRDseedV1::~AliTRDseedV1()
167{
168 //
169 // Destructor. The RecoParam object belongs to the underlying tracker.
170 //
171
29b87567 172 //printf("I-AliTRDseedV1::~AliTRDseedV1() : Owner[%s]\n", IsOwner()?"YES":"NO");
e4f2f73d 173
e3cf3d02 174 if(IsOwner()) {
8d2bec9e 175 for(int itb=0; itb<kNclusters; itb++){
29b87567 176 if(!fClusters[itb]) continue;
177 //AliInfo(Form("deleting c %p @ %d", fClusters[itb], itb));
178 delete fClusters[itb];
4d6aee34 179 fClusters[itb] = NULL;
29b87567 180 }
e3cf3d02 181 }
e4f2f73d 182}
183
184//____________________________________________________________________
185void AliTRDseedV1::Copy(TObject &ref) const
186{
187 //
188 // Copy function
189 //
190
29b87567 191 //AliInfo("");
192 AliTRDseedV1 &target = (AliTRDseedV1 &)ref;
193
4d6aee34 194 target.fkReconstructor = fkReconstructor;
195 target.fClusterIter = NULL;
e3cf3d02 196 target.fExB = fExB;
197 target.fVD = fVD;
198 target.fT0 = fT0;
199 target.fS2PRF = fS2PRF;
200 target.fDiffL = fDiffL;
201 target.fDiffT = fDiffT;
ae4e8b84 202 target.fClusterIdx = 0;
7c3eecb8 203 target.fErrorMsg = fErrorMsg;
3e778975 204 target.fN = fN;
ae4e8b84 205 target.fDet = fDet;
b25a5e9e 206 target.fPt = fPt;
29b87567 207 target.fdX = fdX;
e3cf3d02 208 target.fX0 = fX0;
209 target.fX = fX;
210 target.fY = fY;
211 target.fZ = fZ;
212 target.fS2Y = fS2Y;
213 target.fS2Z = fS2Z;
e3cf3d02 214 target.fChi2 = fChi2;
29b87567 215
8d2bec9e 216 memcpy(target.fIndexes, fIndexes, kNclusters*sizeof(Int_t));
217 memcpy(target.fClusters, fClusters, kNclusters*sizeof(AliTRDcluster*));
2eb10c34 218 memcpy(target.fPad, fPad, 4*sizeof(Float_t));
e3cf3d02 219 target.fYref[0] = fYref[0]; target.fYref[1] = fYref[1];
220 target.fZref[0] = fZref[0]; target.fZref[1] = fZref[1];
221 target.fYfit[0] = fYfit[0]; target.fYfit[1] = fYfit[1];
222 target.fZfit[0] = fZfit[0]; target.fZfit[1] = fZfit[1];
8d2bec9e 223 memcpy(target.fdEdx, fdEdx, kNslices*sizeof(Float_t));
e3cf3d02 224 memcpy(target.fProb, fProb, AliPID::kSPECIES*sizeof(Float_t));
225 memcpy(target.fLabels, fLabels, 3*sizeof(Int_t));
16cca13f 226 memcpy(target.fRefCov, fRefCov, 7*sizeof(Double_t));
68f9b6bd 227 target.fC[0] = fC[0]; target.fC[1] = fC[1];
e3cf3d02 228 memcpy(target.fCov, fCov, 3*sizeof(Double_t));
29b87567 229
e3cf3d02 230 TObject::Copy(ref);
e4f2f73d 231}
232
0906e73e 233
234//____________________________________________________________
4ecadb52 235void AliTRDseedV1::Init(const AliRieman *rieman)
236{
237// Initialize this tracklet using the riemann fit information
238
239
240 fZref[0] = rieman->GetZat(fX0);
241 fZref[1] = rieman->GetDZat(fX0);
242 fYref[0] = rieman->GetYat(fX0);
243 fYref[1] = rieman->GetDYat(fX0);
244 if(fkReconstructor && fkReconstructor->IsHLT()){
245 fRefCov[0] = 1;
246 fRefCov[2] = 10;
247 }else{
248 fRefCov[0] = rieman->GetErrY(fX0);
249 fRefCov[2] = rieman->GetErrZ(fX0);
250 }
251 fC[0] = rieman->GetC();
252 fChi2 = rieman->GetChi2();
253}
254
255
256//____________________________________________________________
f3d3af1b 257Bool_t AliTRDseedV1::Init(AliTRDtrackV1 *track)
0906e73e 258{
259// Initialize this tracklet using the track information
260//
261// Parameters:
262// track - the TRD track used to initialize the tracklet
263//
264// Detailed description
265// The function sets the starting point and direction of the
266// tracklet according to the information from the TRD track.
267//
268// Caution
269// The TRD track has to be propagated to the beginning of the
270// chamber where the tracklet will be constructed
271//
272
29b87567 273 Double_t y, z;
274 if(!track->GetProlongation(fX0, y, z)) return kFALSE;
16cca13f 275 Update(track);
29b87567 276 return kTRUE;
0906e73e 277}
278
bcb6fb78 279
e3cf3d02 280//_____________________________________________________________________________
980d5a2a 281void AliTRDseedV1::Reset(Option_t *opt)
e3cf3d02 282{
980d5a2a 283//
284// Reset seed. If option opt="c" is given only cluster arrays are cleared.
285//
286 for(Int_t ic=kNclusters; ic--;) fIndexes[ic] = -1;
287 memset(fClusters, 0, kNclusters*sizeof(AliTRDcluster*));
560e5c05 288 fN=0; SetBit(kRowCross, kFALSE);
980d5a2a 289 if(strcmp(opt, "c")==0) return;
290
e3cf3d02 291 fExB=0.;fVD=0.;fT0=0.;fS2PRF=0.;
292 fDiffL=0.;fDiffT=0.;
3e778975 293 fClusterIdx=0;
7c3eecb8 294 fErrorMsg = 0;
dd8059a8 295 fDet=-1;
b25a5e9e 296 fPt=0.;
e3cf3d02 297 fdX=0.;fX0=0.; fX=0.; fY=0.; fZ=0.;
298 fS2Y=0.; fS2Z=0.;
68f9b6bd 299 fC[0]=0.; fC[1]=0.;
300 fChi2 = 0.;
e3cf3d02 301
2eb10c34 302 memset(fPad, 0, 4*sizeof(Float_t));
e3cf3d02 303 fYref[0] = 0.; fYref[1] = 0.;
304 fZref[0] = 0.; fZref[1] = 0.;
305 fYfit[0] = 0.; fYfit[1] = 0.;
306 fZfit[0] = 0.; fZfit[1] = 0.;
8d2bec9e 307 memset(fdEdx, 0, kNslices*sizeof(Float_t));
e3cf3d02 308 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = -1.;
309 fLabels[0]=-1; fLabels[1]=-1; // most freq MC labels
310 fLabels[2]=0; // number of different labels for tracklet
16cca13f 311 memset(fRefCov, 0, 7*sizeof(Double_t));
e3cf3d02 312 // covariance matrix [diagonal]
313 // default sy = 200um and sz = 2.3 cm
314 fCov[0] = 4.e-4; fCov[1] = 0.; fCov[2] = 5.3;
315}
316
bcb6fb78 317//____________________________________________________________________
16cca13f 318void AliTRDseedV1::Update(const AliTRDtrackV1 *trk)
b1957d3c 319{
320 // update tracklet reference position from the TRD track
b1957d3c 321
e3cf3d02 322 Double_t fSnp = trk->GetSnp();
323 Double_t fTgl = trk->GetTgl();
b25a5e9e 324 fPt = trk->Pt();
bfd20868 325 Double_t norm =1./TMath::Sqrt((1.-fSnp)*(1.+fSnp));
1fd9389f 326 fYref[1] = fSnp*norm;
327 fZref[1] = fTgl*norm;
b1957d3c 328 SetCovRef(trk->GetCovariance());
329
330 Double_t dx = trk->GetX() - fX0;
331 fYref[0] = trk->GetY() - dx*fYref[1];
332 fZref[0] = trk->GetZ() - dx*fZref[1];
333}
334
e3cf3d02 335//_____________________________________________________________________________
336void AliTRDseedV1::UpdateUsed()
337{
338 //
f29f13a6 339 // Calculate number of used clusers in the tracklet
e3cf3d02 340 //
341
3e778975 342 Int_t nused = 0, nshared = 0;
8d2bec9e 343 for (Int_t i = kNclusters; i--; ) {
e3cf3d02 344 if (!fClusters[i]) continue;
3e778975 345 if(fClusters[i]->IsUsed()){
346 nused++;
347 } else if(fClusters[i]->IsShared()){
348 if(IsStandAlone()) nused++;
349 else nshared++;
350 }
e3cf3d02 351 }
3e778975 352 SetNUsed(nused);
353 SetNShared(nshared);
e3cf3d02 354}
355
356//_____________________________________________________________________________
357void AliTRDseedV1::UseClusters()
358{
359 //
360 // Use clusters
361 //
f29f13a6 362 // In stand alone mode:
363 // Clusters which are marked as used or shared from another track are
364 // removed from the tracklet
365 //
366 // In barrel mode:
367 // - Clusters which are used by another track become shared
368 // - Clusters which are attached to a kink track become shared
369 //
e3cf3d02 370 AliTRDcluster **c = &fClusters[0];
8d2bec9e 371 for (Int_t ic=kNclusters; ic--; c++) {
e3cf3d02 372 if(!(*c)) continue;
f29f13a6 373 if(IsStandAlone()){
374 if((*c)->IsShared() || (*c)->IsUsed()){
b82b4de1 375 if((*c)->IsShared()) SetNShared(GetNShared()-1);
376 else SetNUsed(GetNUsed()-1);
4d6aee34 377 (*c) = NULL;
f29f13a6 378 fIndexes[ic] = -1;
3e778975 379 SetN(GetN()-1);
3e778975 380 continue;
f29f13a6 381 }
3e778975 382 } else {
f29f13a6 383 if((*c)->IsUsed() || IsKink()){
3e778975 384 (*c)->SetShared();
385 continue;
f29f13a6 386 }
387 }
388 (*c)->Use();
e3cf3d02 389 }
390}
391
392
f29f13a6 393
b1957d3c 394//____________________________________________________________________
bcb6fb78 395void AliTRDseedV1::CookdEdx(Int_t nslices)
396{
397// Calculates average dE/dx for all slices and store them in the internal array fdEdx.
398//
399// Parameters:
400// nslices : number of slices for which dE/dx should be calculated
401// Output:
402// store results in the internal array fdEdx. This can be accessed with the method
403// AliTRDseedV1::GetdEdx()
404//
405// Detailed description
406// Calculates average dE/dx for all slices. Depending on the PID methode
407// the number of slices can be 3 (LQ) or 8(NN).
3ee48d6e 408// The calculation of dQ/dl are done using the tracklet fit results (see AliTRDseedV1::GetdQdl(Int_t))
bcb6fb78 409//
410// The following effects are included in the calculation:
411// 1. calibration values for t0 and vdrift (using x coordinate to calculate slice)
412// 2. cluster sharing (optional see AliTRDrecoParam::SetClusterSharing())
413// 3. cluster size
414//
415
8d2bec9e 416 memset(fdEdx, 0, kNslices*sizeof(Float_t));
e73abf77 417 const Double_t kDriftLength = (.5 * AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
29b87567 418
0d80a563 419 AliTRDcluster *c(NULL);
420 for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++){
421 if(!(c = fClusters[ic]) && !(c = fClusters[ic+kNtb])) continue;
422 Float_t dx = TMath::Abs(fX0 - c->GetX());
423
29b87567 424 // Filter clusters for dE/dx calculation
0d80a563 425
29b87567 426 // 1.consider calibration effects for slice determination
0d80a563 427 Int_t slice;
428 if(dx<kDriftLength){ // TODO should be replaced by c->IsInChamber()
429 slice = Int_t(dx * nslices / kDriftLength);
430 } else slice = c->GetX() < fX0 ? nslices-1 : 0;
431
e73abf77 432
29b87567 433 // 2. take sharing into account
3e778975 434 Float_t w = /*c->IsShared() ? .5 :*/ 1.;
0d80a563 435
29b87567 436 // 3. take into account large clusters TODO
437 //w *= c->GetNPads() > 3 ? .8 : 1.;
85b917f6 438
0d80a563 439 //CHECK !!!
440 fdEdx[slice] += w * GetdQdl(ic); //fdQdl[ic];
441 } // End of loop over clusters
bcb6fb78 442}
443
e3cf3d02 444//_____________________________________________________________________________
445void AliTRDseedV1::CookLabels()
446{
447 //
448 // Cook 2 labels for seed
449 //
450
451 Int_t labels[200];
452 Int_t out[200];
453 Int_t nlab = 0;
8d2bec9e 454 for (Int_t i = 0; i < kNclusters; i++) {
e3cf3d02 455 if (!fClusters[i]) continue;
456 for (Int_t ilab = 0; ilab < 3; ilab++) {
457 if (fClusters[i]->GetLabel(ilab) >= 0) {
458 labels[nlab] = fClusters[i]->GetLabel(ilab);
459 nlab++;
460 }
461 }
462 }
463
fac58f00 464 fLabels[2] = AliMathBase::Freq(nlab,labels,out,kTRUE);
e3cf3d02 465 fLabels[0] = out[0];
466 if ((fLabels[2] > 1) && (out[3] > 1)) fLabels[1] = out[2];
467}
468
2eb10c34 469//____________________________________________________________
470Float_t AliTRDseedV1::GetAnodeWireOffset(Float_t zt)
471{
4ecadb52 472// Find position inside the amplification cell for reading drift velocity map
473
2eb10c34 474 Float_t d = fPad[3] - zt;
475 if(d<0.){
476 AliError(Form("Fail AnodeWireOffset calculation z0[%+7.2f] zt[%+7.2f] d[%+7.2f].", fPad[3], zt, d));
477 return 0.125;
478 }
479 d -= ((Int_t)(2 * d)) / 2.0;
480 if(d > 0.25) d = 0.5 - d;
481 return d;
482}
483
e3cf3d02 484
d937ad7a 485//____________________________________________________________________
85b917f6 486Float_t AliTRDseedV1::GetCharge(Bool_t useOutliers)
487{
488// Computes total charge attached to tracklet. If "useOutliers" is set clusters
489// which are not in chamber are also used (default false)
490
491 AliTRDcluster *c(NULL); Float_t qt(0.);
492 for(int ic=0; ic<kNclusters; ic++){
493 if(!(c=fClusters[ic])) continue;
494 if(c->IsInChamber() && !useOutliers) continue;
495 qt += TMath::Abs(c->GetQ());
496 }
497 return qt;
498}
499
500//____________________________________________________________________
e1bcf0af 501Bool_t AliTRDseedV1::GetEstimatedCrossPoint(Float_t &x, Float_t &z) const
502{
503// Algorithm to estimate cross point in the x-z plane for pad row cross tracklets.
504// Returns true in case of success.
505 if(!IsRowCross()) return kFALSE;
506
507 x=0.; z=0.;
508 AliTRDcluster *c(NULL);
509 // Find radial range for first row
510 Float_t x1[] = {0., 1.e3};
511 for(int ic=0; ic<kNtb; ic++){
512 if(!(c=fClusters[ic])) continue;
513 if(!c->IsInChamber()) continue;
514 if(c->GetX() <= x1[1]) x1[1] = c->GetX();
515 if(c->GetX() >= x1[0]) x1[0] = c->GetX();
516 z=c->GetZ();
517 }
518 if((x1[0] - x1[1])<1.e-5) return kFALSE;
519
520 // Find radial range for second row
521 Bool_t kZ(kFALSE);
522 Float_t x2[] = {0., 1.e3};
523 for(int ic=kNtb; ic<kNclusters; ic++){
524 if(!(c=fClusters[ic])) continue;
525 if(!c->IsInChamber()) continue;
526 if(c->GetX() <= x2[1]) x2[1] = c->GetX();
527 if(c->GetX() >= x2[0]) x2[0] = c->GetX();
528 if(!kZ){
529 z+=c->GetZ();
530 z*=0.5;
531 kZ=kTRUE;
532 }
533 }
534 if((x2[0] - x2[1])<1.e-5) return kFALSE;
535
536 // Find intersection of the 2 radial regions
537 x = 0.5*((x1[0]+x1[1] > x2[0]+x2[1]) ? (x1[1]+x2[0]) : (x1[0]+x2[1]));
538 return kTRUE;
539}
540
541//____________________________________________________________________
0b433f72 542Float_t AliTRDseedV1::GetdQdl(Int_t ic, Float_t *dl) const
bcb6fb78 543{
3ee48d6e 544// Using the linear approximation of the track inside one TRD chamber (TRD tracklet)
545// the charge per unit length can be written as:
546// BEGIN_LATEX
500851ab 547// #frac{dq}{dl} = #frac{q_{c}}{dx * #sqrt{1 + #(){#frac{dy}{dx}}^{2}_{fit} + #(){#frac{dz}{dx}}^{2}_{ref}}}
3ee48d6e 548// END_LATEX
549// where qc is the total charge collected in the current time bin and dx is the length
0b433f72 550// of the time bin.
551// The following correction are applied :
552// - charge : pad row cross corrections
553// [diffusion and TRF assymetry] TODO
554// - dx : anisochronity, track inclination - see Fit and AliTRDcluster::GetXloc()
555// and AliTRDcluster::GetYloc() for the effects taken into account
3ee48d6e 556//
0fa1a8ee 557//Begin_Html
558//<img src="TRD/trackletDQDT.gif">
559//End_Html
560// In the picture the energy loss measured on the tracklet as a function of drift time [left] and respectively
561// drift length [right] for different particle species is displayed.
3ee48d6e 562// Author : Alex Bercuci <A.Bercuci@gsi.de>
563//
564 Float_t dq = 0.;
5d401b45 565 // check whether both clusters are inside the chamber
566 Bool_t hasClusterInChamber = kFALSE;
567 if(fClusters[ic] && fClusters[ic]->IsInChamber()){
568 hasClusterInChamber = kTRUE;
1742f24c 569 dq += TMath::Abs(fClusters[ic]->GetQ());
b30d8c09 570 }
571 if(fClusters[ic+kNtb] && fClusters[ic+kNtb]->IsInChamber()){
5d401b45 572 hasClusterInChamber = kTRUE;
573 dq += TMath::Abs(fClusters[ic+kNtb]->GetQ());
1742f24c 574 }
5d401b45 575 if(!hasClusterInChamber) return 0.;
0b433f72 576 if(dq<1.e-3) return 0.;
3ee48d6e 577
a2abcbc5 578 Double_t dx = fdX;
579 if(ic-1>=0 && ic+1<kNtb){
580 Float_t x2(0.), x1(0.);
5d401b45 581 // try to estimate upper radial position (find the cluster which is inside the chamber)
582 if(fClusters[ic-1] && fClusters[ic-1]->IsInChamber()) x2 = fClusters[ic-1]->GetX();
583 else if(fClusters[ic-1+kNtb] && fClusters[ic-1+kNtb]->IsInChamber()) x2 = fClusters[ic-1+kNtb]->GetX();
584 else if(fClusters[ic] && fClusters[ic]->IsInChamber()) x2 = fClusters[ic]->GetX()+fdX;
a2abcbc5 585 else x2 = fClusters[ic+kNtb]->GetX()+fdX;
5d401b45 586 // try to estimate lower radial position (find the cluster which is inside the chamber)
587 if(fClusters[ic+1] && fClusters[ic+1]->IsInChamber()) x1 = fClusters[ic+1]->GetX();
588 else if(fClusters[ic+1+kNtb] && fClusters[ic+1+kNtb]->IsInChamber()) x1 = fClusters[ic+1+kNtb]->GetX();
589 else if(fClusters[ic] && fClusters[ic]->IsInChamber()) x1 = fClusters[ic]->GetX()-fdX;
a2abcbc5 590 else x1 = fClusters[ic+kNtb]->GetX()-fdX;
591
592 dx = .5*(x2 - x1);
593 }
0b433f72 594 dx *= TMath::Sqrt(1. + fYfit[1]*fYfit[1] + fZref[1]*fZref[1]);
0b433f72 595 if(dl) (*dl) = dx;
283604d2 596 if(dx>1.e-9) return dq/dx;
597 else return 0.;
bcb6fb78 598}
599
0b433f72 600//____________________________________________________________
601Float_t AliTRDseedV1::GetMomentum(Float_t *err) const
602{
603// Returns momentum of the track after update with the current tracklet as:
604// BEGIN_LATEX
605// p=#frac{1}{1/p_{t}} #sqrt{1+tgl^{2}}
606// END_LATEX
607// and optionally the momentum error (if err is not null).
608// The estimated variance of the momentum is given by:
609// BEGIN_LATEX
610// #sigma_{p}^{2} = (#frac{dp}{dp_{t}})^{2} #sigma_{p_{t}}^{2}+(#frac{dp}{dtgl})^{2} #sigma_{tgl}^{2}+2#frac{dp}{dp_{t}}#frac{dp}{dtgl} cov(tgl,1/p_{t})
611// END_LATEX
612// which can be simplified to
613// BEGIN_LATEX
614// #sigma_{p}^{2} = p^{2}p_{t}^{4}tgl^{2}#sigma_{tgl}^{2}-2p^{2}p_{t}^{3}tgl cov(tgl,1/p_{t})+p^{2}p_{t}^{2}#sigma_{1/p_{t}}^{2}
615// END_LATEX
616//
617
618 Double_t p = fPt*TMath::Sqrt(1.+fZref[1]*fZref[1]);
619 Double_t p2 = p*p;
620 Double_t tgl2 = fZref[1]*fZref[1];
621 Double_t pt2 = fPt*fPt;
622 if(err){
623 Double_t s2 =
624 p2*tgl2*pt2*pt2*fRefCov[4]
625 -2.*p2*fZref[1]*fPt*pt2*fRefCov[5]
626 +p2*pt2*fRefCov[6];
627 (*err) = TMath::Sqrt(s2);
628 }
629 return p;
630}
631
b453ef55 632//____________________________________________________________________
633Float_t AliTRDseedV1::GetOccupancyTB() const
634{
635// Returns procentage of TB occupied by clusters
636
637 Int_t n(0);
638 AliTRDcluster *c(NULL);
639 for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++){
640 if(!(c = fClusters[ic]) && !(c = fClusters[ic+kNtb])) continue;
641 n++;
642 }
643
644 return Float_t(n)/AliTRDtrackerV1::GetNTimeBins();
645}
0b433f72 646
0906e73e 647//____________________________________________________________________
3e778975 648Float_t* AliTRDseedV1::GetProbability(Bool_t force)
0906e73e 649{
3e778975 650 if(!force) return &fProb[0];
4d6aee34 651 if(!CookPID()) return NULL;
3e778975 652 return &fProb[0];
653}
654
655//____________________________________________________________
656Bool_t AliTRDseedV1::CookPID()
657{
0906e73e 658// Fill probability array for tracklet from the DB.
659//
660// Parameters
661//
662// Output
4d6aee34 663// returns pointer to the probability array and NULL if missing DB access
0906e73e 664//
2a3191bb 665// Retrieve PID probabilities for e+-, mu+-, K+-, pi+- and p+- from the DB according to tracklet information:
666// - estimated momentum at tracklet reference point
667// - dE/dx measurements
668// - tracklet length
669// - TRD layer
670// According to the steering settings specified in the reconstruction one of the following methods are used
671// - Neural Network [default] - option "nn"
672// - 2D Likelihood - option "!nn"
0906e73e 673
0906e73e 674 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
675 if (!calibration) {
676 AliError("No access to calibration data");
3e778975 677 return kFALSE;
0906e73e 678 }
679
4d6aee34 680 if (!fkReconstructor) {
3a039a31 681 AliError("Reconstructor not set.");
3e778975 682 return kFALSE;
4ba1d6ae 683 }
684
0906e73e 685 // Retrieve the CDB container class with the parametric detector response
4d6aee34 686 const AliTRDCalPID *pd = calibration->GetPIDObject(fkReconstructor->GetPIDMethod());
0906e73e 687 if (!pd) {
688 AliError("No access to AliTRDCalPID object");
3e778975 689 return kFALSE;
0906e73e 690 }
10f75631 691
29b87567 692 // calculate tracklet length TO DO
560e5c05 693 Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())/ TMath::Sqrt((1.0 - GetSnp()*GetSnp()) / (1.0 + GetTgl()*GetTgl()));
0906e73e 694
695 //calculate dE/dx
9ded305e 696 CookdEdx(AliTRDCalPID::kNSlicesNN);
697 AliDebug(4, Form("p=%6.4f[GeV/c] dEdx{%7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f} l=%4.2f[cm]", GetMomentum(), fdEdx[0], fdEdx[1], fdEdx[2], fdEdx[3], fdEdx[4], fdEdx[5], fdEdx[6], fdEdx[7], length));
0217fcd0 698
0906e73e 699 // Sets the a priori probabilities
11d80e40 700 Bool_t kPIDNN(fkReconstructor->GetPIDMethod()==AliTRDpidUtil::kNN);
f83cd814 701 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++)
11d80e40 702 fProb[ispec] = pd->GetProbability(ispec, GetMomentum(), &fdEdx[0], length, kPIDNN?GetPlane():fkReconstructor->GetRecoParam()->GetPIDLQslices());
f301a656 703
3e778975 704 return kTRUE;
0906e73e 705}
706
707//____________________________________________________________________
e4f2f73d 708Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
709{
710 //
711 // Returns a quality measurement of the current seed
712 //
713
dd8059a8 714 Float_t zcorr = kZcorr ? GetTilt() * (fZfit[0] - fZref[0]) : 0.;
29b87567 715 return
3e778975 716 .5 * TMath::Abs(18.0 - GetN())
29b87567 717 + 10.* TMath::Abs(fYfit[1] - fYref[1])
718 + 5. * TMath::Abs(fYfit[0] - fYref[0] + zcorr)
dd8059a8 719 + 2. * TMath::Abs(fZfit[0] - fZref[0]) / GetPadLength();
e4f2f73d 720}
721
722//____________________________________________________________________
d937ad7a 723void AliTRDseedV1::GetCovAt(Double_t x, Double_t *cov) const
0906e73e 724{
d937ad7a 725// Computes covariance in the y-z plane at radial point x (in tracking coordinates)
726// and returns the results in the preallocated array cov[3] as :
727// cov[0] = Var(y)
728// cov[1] = Cov(yz)
729// cov[2] = Var(z)
730//
731// Details
732//
733// For the linear transformation
734// BEGIN_LATEX
735// Y = T_{x} X^{T}
736// END_LATEX
737// The error propagation has the general form
738// BEGIN_LATEX
739// C_{Y} = T_{x} C_{X} T_{x}^{T}
740// END_LATEX
741// We apply this formula 2 times. First to calculate the covariance of the tracklet
742// at point x we consider:
743// BEGIN_LATEX
744// T_{x} = (1 x); X=(y0 dy/dx); C_{X}=#(){#splitline{Var(y0) Cov(y0, dy/dx)}{Cov(y0, dy/dx) Var(dy/dx)}}
745// END_LATEX
746// and secondly to take into account the tilt angle
747// BEGIN_LATEX
748// T_{#alpha} = #(){#splitline{cos(#alpha) __ sin(#alpha)}{-sin(#alpha) __ cos(#alpha)}}; X=(y z); C_{X}=#(){#splitline{Var(y) 0}{0 Var(z)}}
749// END_LATEX
750//
751// using simple trigonometrics one can write for this last case
752// BEGIN_LATEX
753// 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})}}
754// END_LATEX
755// which can be aproximated for small alphas (2 deg) with
756// BEGIN_LATEX
757// 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}}}
758// END_LATEX
759//
760// before applying the tilt rotation we also apply systematic uncertainties to the tracklet
761// position which can be tunned from outside via the AliTRDrecoParam::SetSysCovMatrix(). They might
762// account for extra misalignment/miscalibration uncertainties.
763//
764// Author :
765// Alex Bercuci <A.Bercuci@gsi.de>
766// Date : Jan 8th 2009
767//
b1957d3c 768
769
d937ad7a 770 Double_t xr = fX0-x;
771 Double_t sy2 = fCov[0] +2.*xr*fCov[1] + xr*xr*fCov[2];
b72f4eaf 772 Double_t sz2 = fS2Z;
773 //GetPadLength()*GetPadLength()/12.;
0906e73e 774
d937ad7a 775 // insert systematic uncertainties
4d6aee34 776 if(fkReconstructor){
bb2db46c 777 Double_t sys[15]; memset(sys, 0, 15*sizeof(Double_t));
4d6aee34 778 fkReconstructor->GetRecoParam()->GetSysCovMatrix(sys);
bb2db46c 779 sy2 += sys[0];
780 sz2 += sys[1];
781 }
2eb10c34 782
783 // rotate covariance matrix if no RC
784 if(!IsRowCross()){
785 Double_t t2 = GetTilt()*GetTilt();
786 Double_t correction = 1./(1. + t2);
787 cov[0] = (sy2+t2*sz2)*correction;
788 cov[1] = GetTilt()*(sz2 - sy2)*correction;
789 cov[2] = (t2*sy2+sz2)*correction;
790 } else {
791 cov[0] = sy2; cov[1] = 0.; cov[2] = sy2;
792 }
793
794 AliDebug(4, Form("C(%6.1f %+6.3f %6.1f) RC[%c]", 1.e4*TMath::Sqrt(cov[0]), cov[1], 1.e4*TMath::Sqrt(cov[2]), IsRowCross()?'y':'n'));
d937ad7a 795}
eb38ed55 796
bb2db46c 797//____________________________________________________________
66765e8e 798Int_t AliTRDseedV1::GetCovSqrt(const Double_t * const c, Double_t *d)
bb2db46c 799{
800// Helper function to calculate the square root of the covariance matrix.
801// The input matrix is stored in the vector c and the result in the vector d.
41b7c7b6 802// Both arrays have to be initialized by the user with at least 3 elements. Return negative in case of failure.
bb2db46c 803//
ec3f0161 804// For calculating the square root of the symmetric matrix c
805// the following relation is used:
bb2db46c 806// BEGIN_LATEX
ec3f0161 807// C^{1/2} = VD^{1/2}V^{-1}
bb2db46c 808// END_LATEX
41b7c7b6 809// with V being the matrix with the n eigenvectors as columns.
ec3f0161 810// In case C is symmetric the followings are true:
811// - matrix D is diagonal with the diagonal given by the eigenvalues of C
41b7c7b6 812// - V = V^{-1}
bb2db46c 813//
814// Author A.Bercuci <A.Bercuci@gsi.de>
815// Date Mar 19 2009
816
66765e8e 817 const Double_t kZero(1.e-20);
4d6aee34 818 Double_t l[2], // eigenvalues
819 v[3]; // eigenvectors
bb2db46c 820 // the secular equation and its solution :
821 // (c[0]-L)(c[2]-L)-c[1]^2 = 0
822 // L^2 - L*Tr(c)+DET(c) = 0
823 // L12 = [Tr(c) +- sqrt(Tr(c)^2-4*DET(c))]/2
4d6aee34 824 Double_t tr = c[0]+c[2], // trace
825 det = c[0]*c[2]-c[1]*c[1]; // determinant
66765e8e 826 if(TMath::Abs(det)<kZero) return 1;
4d6aee34 827 Double_t dd = TMath::Sqrt(tr*tr - 4*det);
66765e8e 828 l[0] = .5*(tr + dd*(c[0]>c[2]?-1.:1.));
829 l[1] = .5*(tr + dd*(c[0]>c[2]?1.:-1.));
830 if(l[0]<kZero || l[1]<kZero) return 2;
41b7c7b6 831 // the sym V matrix
832 // | v00 v10|
833 // | v10 v11|
66765e8e 834 Double_t den = (l[0]-c[0])*(l[0]-c[0])+c[1]*c[1];
835 if(den<kZero){ // almost diagonal
836 v[0] = TMath::Sign(0., c[1]);
837 v[1] = TMath::Sign(1., (l[0]-c[0]));
838 v[2] = TMath::Sign(0., c[1]*(l[0]-c[0])*(l[1]-c[2]));
839 } else {
840 Double_t tmp = 1./TMath::Sqrt(den);
841 v[0] = c[1]* tmp;
842 v[1] = (l[0]-c[0])*tmp;
843 if(TMath::Abs(l[1]-c[2])<kZero) v[2] = TMath::Sign(v[0]*(l[0]-c[0])/kZero, (l[1]-c[2]));
844 else v[2] = v[0]*(l[0]-c[0])/(l[1]-c[2]);
845 }
41b7c7b6 846 // the VD^{1/2}V is:
4d6aee34 847 l[0] = TMath::Sqrt(l[0]); l[1] = TMath::Sqrt(l[1]);
848 d[0] = v[0]*v[0]*l[0]+v[1]*v[1]*l[1];
849 d[1] = v[0]*v[1]*l[0]+v[1]*v[2]*l[1];
850 d[2] = v[1]*v[1]*l[0]+v[2]*v[2]*l[1];
bb2db46c 851
66765e8e 852 return 0;
bb2db46c 853}
854
855//____________________________________________________________
4d6aee34 856Double_t AliTRDseedV1::GetCovInv(const Double_t * const c, Double_t *d)
bb2db46c 857{
858// Helper function to calculate the inverse of the covariance matrix.
859// The input matrix is stored in the vector c and the result in the vector d.
860// Both arrays have to be initialized by the user with at least 3 elements
861// The return value is the determinant or 0 in case of singularity.
862//
863// Author A.Bercuci <A.Bercuci@gsi.de>
864// Date Mar 19 2009
865
4d6aee34 866 Double_t det = c[0]*c[2] - c[1]*c[1];
867 if(TMath::Abs(det)<1.e-20) return 0.;
868 Double_t invDet = 1./det;
869 d[0] = c[2]*invDet;
870 d[1] =-c[1]*invDet;
871 d[2] = c[0]*invDet;
872 return det;
bb2db46c 873}
0906e73e 874
d937ad7a 875//____________________________________________________________________
b72f4eaf 876UShort_t AliTRDseedV1::GetVolumeId() const
877{
4ecadb52 878// Returns geometry volume id by delegation
879
fbe11be7 880 for(Int_t ic(0);ic<kNclusters; ic++){
881 if(fClusters[ic]) return fClusters[ic]->GetVolumeId();
882 }
883 return 0;
b72f4eaf 884}
885
886
887//____________________________________________________________________
e3cf3d02 888void AliTRDseedV1::Calibrate()
d937ad7a 889{
e3cf3d02 890// Retrieve calibration and position parameters from OCDB.
891// The following information are used
d937ad7a 892// - detector index
e3cf3d02 893// - column and row position of first attached cluster. If no clusters are attached
894// to the tracklet a random central chamber position (c=70, r=7) will be used.
895//
896// The following information is cached in the tracklet
897// t0 (trigger delay)
898// drift velocity
899// PRF width
900// omega*tau = tg(a_L)
901// diffusion coefficients (longitudinal and transversal)
d937ad7a 902//
903// Author :
904// Alex Bercuci <A.Bercuci@gsi.de>
905// Date : Jan 8th 2009
906//
eb38ed55 907
d937ad7a 908 AliCDBManager *cdb = AliCDBManager::Instance();
909 if(cdb->GetRun() < 0){
910 AliError("OCDB manager not properly initialized");
911 return;
912 }
0906e73e 913
e3cf3d02 914 AliTRDcalibDB *calib = AliTRDcalibDB::Instance();
915 AliTRDCalROC *vdROC = calib->GetVdriftROC(fDet),
916 *t0ROC = calib->GetT0ROC(fDet);;
917 const AliTRDCalDet *vdDet = calib->GetVdriftDet();
918 const AliTRDCalDet *t0Det = calib->GetT0Det();
d937ad7a 919
920 Int_t col = 70, row = 7;
921 AliTRDcluster **c = &fClusters[0];
3e778975 922 if(GetN()){
d937ad7a 923 Int_t ic = 0;
8d2bec9e 924 while (ic<kNclusters && !(*c)){ic++; c++;}
d937ad7a 925 if(*c){
926 col = (*c)->GetPadCol();
927 row = (*c)->GetPadRow();
928 }
929 }
3a039a31 930
e17f4785 931 fT0 = (t0Det->GetValue(fDet) + t0ROC->GetValue(col,row)) / AliTRDCommonParam::Instance()->GetSamplingFrequency();
e3cf3d02 932 fVD = vdDet->GetValue(fDet) * vdROC->GetValue(col, row);
933 fS2PRF = calib->GetPRFWidth(fDet, col, row); fS2PRF *= fS2PRF;
934 fExB = AliTRDCommonParam::Instance()->GetOmegaTau(fVD);
935 AliTRDCommonParam::Instance()->GetDiffCoeff(fDiffL,
936 fDiffT, fVD);
903326c1 937 AliDebug(4, Form("Calibration params for Det[%3d] Col[%3d] Row[%2d]\n t0[%f] vd[%f] s2PRF[%f] ExB[%f] Dl[%f] Dt[%f]", fDet, col, row, fT0, fVD, fS2PRF, fExB, fDiffL, fDiffT));
938
939
e3cf3d02 940 SetBit(kCalib, kTRUE);
0906e73e 941}
942
0906e73e 943//____________________________________________________________________
29b87567 944void AliTRDseedV1::SetOwner()
0906e73e 945{
29b87567 946 //AliInfo(Form("own [%s] fOwner[%s]", own?"YES":"NO", fOwner?"YES":"NO"));
947
948 if(TestBit(kOwner)) return;
8d2bec9e 949 for(int ic=0; ic<kNclusters; ic++){
29b87567 950 if(!fClusters[ic]) continue;
951 fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
952 }
953 SetBit(kOwner);
0906e73e 954}
955
eb2b4f91 956//____________________________________________________________
4ecadb52 957void AliTRDseedV1::SetPadPlane(AliTRDpadPlane * const p)
eb2b4f91 958{
959// Shortcut method to initialize pad geometry.
2eb10c34 960 fPad[0] = p->GetLengthIPad();
961 fPad[1] = p->GetWidthIPad();
962 fPad[2] = TMath::Tan(TMath::DegToRad()*p->GetTiltingAngle());
963 fPad[3] = p->GetRow0() + p->GetAnodeWireOffset();
eb2b4f91 964}
965
966
e4f2f73d 967//____________________________________________________________________
4d6aee34 968Bool_t AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *const chamber, Bool_t tilt)
e4f2f73d 969{
1fd9389f 970//
971// Projective algorithm to attach clusters to seeding tracklets. The following steps are performed :
972// 1. Collapse x coordinate for the full detector plane
973// 2. truncated mean on y (r-phi) direction
974// 3. purge clusters
975// 4. truncated mean on z direction
976// 5. purge clusters
977//
978// Parameters
979// - chamber : pointer to tracking chamber container used to search the tracklet
980// - tilt : switch for tilt correction during road building [default true]
981// Output
982// - true : if tracklet found successfully. Failure can happend because of the following:
983// -
984// Detailed description
985//
986// We start up by defining the track direction in the xy plane and roads. The roads are calculated based
8a7ff53c 987// on tracking information (variance in the r-phi direction) and estimated variance of the standard
988// clusters (see AliTRDcluster::SetSigmaY2()) corrected for tilt (see GetCovAt()). From this the road is
989// BEGIN_LATEX
500851ab 990// r_{y} = 3*#sqrt{12*(#sigma^{2}_{Trk}(y) + #frac{#sigma^{2}_{cl}(y) + tg^{2}(#alpha_{L})#sigma^{2}_{cl}(z)}{1+tg^{2}(#alpha_{L})})}
8a7ff53c 991// r_{z} = 1.5*L_{pad}
992// END_LATEX
1fd9389f 993//
4b755889 994// Author : Alexandru Bercuci <A.Bercuci@gsi.de>
995// Debug : level >3
1fd9389f 996
fc0882f3 997 const AliTRDrecoParam* const recoParam = fkReconstructor->GetRecoParam(); //the dynamic cast in GetRecoParam is slow, so caching the pointer to it
998
999 if(!recoParam){
560e5c05 1000 AliError("Tracklets can not be used without a valid RecoParam.");
29b87567 1001 return kFALSE;
1002 }
b1957d3c 1003 // Initialize reco params for this tracklet
1004 // 1. first time bin in the drift region
a2abcbc5 1005 Int_t t0 = 14;
fc0882f3 1006 Int_t kClmin = Int_t(recoParam->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
29b87567 1007
fc0882f3 1008 Double_t sysCov[5]; recoParam->GetSysCovMatrix(sysCov);
8a7ff53c 1009 Double_t s2yTrk= fRefCov[0],
1010 s2yCl = 0.,
1011 s2zCl = GetPadLength()*GetPadLength()/12.,
1012 syRef = TMath::Sqrt(s2yTrk),
1013 t2 = GetTilt()*GetTilt();
29b87567 1014 //define roads
fc0882f3 1015 Double_t kroady = 1., //recoParam->GetRoad1y();
1016 kroadz = GetPadLength() * recoParam->GetRoadzMultiplicator() + 1.;
8a7ff53c 1017 // define probing cluster (the perfect cluster) and default calibration
1018 Short_t sig[] = {0, 0, 10, 30, 10, 0,0};
1019 AliTRDcluster cp(fDet, 6, 75, 0, sig, 0);
560e5c05 1020 if(fkReconstructor->IsHLT()) cp.SetRPhiMethod(AliTRDcluster::kCOG);
1021 if(!IsCalibrated()) Calibrate();
8a7ff53c 1022
ee8fb199 1023 AliDebug(4, "");
1024 AliDebug(4, Form("syKalman[%f] rY[%f] rZ[%f]", syRef, kroady, kroadz));
29b87567 1025
1026 // working variables
b1957d3c 1027 const Int_t kNrows = 16;
4b755889 1028 const Int_t kNcls = 3*kNclusters; // buffer size
1029 AliTRDcluster *clst[kNrows][kNcls];
3044dfe5 1030 Bool_t blst[kNrows][kNcls];
4b755889 1031 Double_t cond[4], dx, dy, yt, zt, yres[kNrows][kNcls];
1032 Int_t idxs[kNrows][kNcls], ncl[kNrows], ncls = 0;
b1957d3c 1033 memset(ncl, 0, kNrows*sizeof(Int_t));
4b755889 1034 memset(yres, 0, kNrows*kNcls*sizeof(Double_t));
3044dfe5 1035 memset(blst, 0, kNrows*kNcls*sizeof(Bool_t)); //this is 8 times faster to memset than "memset(clst, 0, kNrows*kNcls*sizeof(AliTRDcluster*))"
b1957d3c 1036
29b87567 1037 // Do cluster projection
4d6aee34 1038 AliTRDcluster *c = NULL;
1039 AliTRDchamberTimeBin *layer = NULL;
b1957d3c 1040 Bool_t kBUFFER = kFALSE;
4b755889 1041 for (Int_t it = 0; it < kNtb; it++) {
b1957d3c 1042 if(!(layer = chamber->GetTB(it))) continue;
29b87567 1043 if(!Int_t(*layer)) continue;
8a7ff53c 1044 // get track projection at layers position
b1957d3c 1045 dx = fX0 - layer->GetX();
1046 yt = fYref[0] - fYref[1] * dx;
1047 zt = fZref[0] - fZref[1] * dx;
8a7ff53c 1048 // get standard cluster error corrected for tilt
1049 cp.SetLocalTimeBin(it);
1050 cp.SetSigmaY2(0.02, fDiffT, fExB, dx, -1./*zt*/, fYref[1]);
d956a643 1051 s2yCl = (cp.GetSigmaY2() + sysCov[0] + t2*s2zCl)/(1.+t2);
8a7ff53c 1052 // get estimated road
1053 kroady = 3.*TMath::Sqrt(12.*(s2yTrk + s2yCl));
1054
ee8fb199 1055 AliDebug(5, Form(" %2d x[%f] yt[%f] zt[%f]", it, dx, yt, zt));
1056
1057 AliDebug(5, Form(" syTrk[um]=%6.2f syCl[um]=%6.2f syClTlt[um]=%6.2f Ry[mm]=%f", 1.e4*TMath::Sqrt(s2yTrk), 1.e4*TMath::Sqrt(cp.GetSigmaY2()), 1.e4*TMath::Sqrt(s2yCl), 1.e1*kroady));
b1957d3c 1058
8a7ff53c 1059 // select clusters
b1957d3c 1060 cond[0] = yt; cond[2] = kroady;
1061 cond[1] = zt; cond[3] = kroadz;
1062 Int_t n=0, idx[6];
1063 layer->GetClusters(cond, idx, n, 6);
1064 for(Int_t ic = n; ic--;){
1065 c = (*layer)[idx[ic]];
1066 dy = yt - c->GetY();
dd8059a8 1067 dy += tilt ? GetTilt() * (c->GetZ() - zt) : 0.;
b1957d3c 1068 // select clusters on a 3 sigmaKalman level
1069/* if(tilt && TMath::Abs(dy) > 3.*syRef){
1070 printf("too large !!!\n");
1071 continue;
1072 }*/
1073 Int_t r = c->GetPadRow();
ee8fb199 1074 AliDebug(5, Form(" -> dy[%f] yc[%f] r[%d]", TMath::Abs(dy), c->GetY(), r));
b1957d3c 1075 clst[r][ncl[r]] = c;
3044dfe5 1076 blst[r][ncl[r]] = kTRUE;
b1957d3c 1077 idxs[r][ncl[r]] = idx[ic];
1078 yres[r][ncl[r]] = dy;
1079 ncl[r]++; ncls++;
1080
4b755889 1081 if(ncl[r] >= kNcls) {
560e5c05 1082 AliWarning(Form("Cluster candidates row[%d] reached buffer limit[%d]. Some may be lost.", r, kNcls));
b1957d3c 1083 kBUFFER = kTRUE;
29b87567 1084 break;
1085 }
1086 }
b1957d3c 1087 if(kBUFFER) break;
29b87567 1088 }
ee8fb199 1089 AliDebug(4, Form("Found %d clusters. Processing ...", ncls));
1090 if(ncls<kClmin){
560e5c05 1091 AliDebug(1, Form("CLUSTERS FOUND %d LESS THAN THRESHOLD %d.", ncls, kClmin));
7c3eecb8 1092 SetErrorMsg(kAttachClFound);
ee8fb199 1093 return kFALSE;
1094 }
1095
b1957d3c 1096 // analyze each row individualy
560e5c05 1097 Bool_t kRowSelection(kFALSE);
1098 Double_t mean[]={1.e3, 1.e3, 1.3}, syDis[]={1.e3, 1.e3, 1.3};
1099 Int_t nrow[] = {0, 0, 0}, rowId[] = {-1, -1, -1}, nr = 0, lr=-1;
1100 TVectorD vdy[3];
1101 for(Int_t ir=0; ir<kNrows; ir++){
b1957d3c 1102 if(!(ncl[ir])) continue;
560e5c05 1103 if(lr>0 && ir-lr != 1){
1104 AliDebug(2, "Rows attached not continuous. Turn on selection.");
1105 kRowSelection=kTRUE;
29b87567 1106 }
560e5c05 1107
ee8fb199 1108 AliDebug(5, Form(" r[%d] n[%d]", ir, ncl[ir]));
b1957d3c 1109 // Evaluate truncated mean on the y direction
560e5c05 1110 if(ncl[ir] < 4) continue;
1111 AliMathBase::EvaluateUni(ncl[ir], yres[ir], mean[nr], syDis[nr], Int_t(ncl[ir]*.8));
4b755889 1112
b1957d3c 1113 // TODO check mean and sigma agains cluster resolution !!
560e5c05 1114 AliDebug(4, Form(" m_%d[%+5.3f (%5.3fs)] s[%f]", nr, mean[nr], TMath::Abs(mean[nr]/syDis[nr]), syDis[nr]));
1115 // remove outliers based on a 3 sigmaDistr level
b1957d3c 1116 Bool_t kFOUND = kFALSE;
1117 for(Int_t ic = ncl[ir]; ic--;){
560e5c05 1118 if(yres[ir][ic] - mean[nr] > 3. * syDis[nr]){
3044dfe5 1119 blst[ir][ic] = kFALSE; continue;
b1957d3c 1120 }
560e5c05 1121 nrow[nr]++; rowId[nr]=ir; kFOUND = kTRUE;
1122 }
1123 if(kFOUND){
1124 vdy[nr].Use(nrow[nr], yres[ir]);
1125 nr++;
b1957d3c 1126 }
b1957d3c 1127 lr = ir; if(nr>=3) break;
29b87567 1128 }
fc0882f3 1129 if(recoParam->GetStreamLevel(AliTRDrecoParam::kTracker) > 3 && fkReconstructor->IsDebugStreaming()){
560e5c05 1130 TTreeSRedirector &cstreamer = *fkReconstructor->GetDebugStream(AliTRDrecoParam::kTracker);
1131 UChar_t stat(0);
1132 if(IsKink()) SETBIT(stat, 1);
1133 if(IsStandAlone()) SETBIT(stat, 2);
1134 cstreamer << "AttachClusters"
1135 << "stat=" << stat
1136 << "det=" << fDet
1137 << "pt=" << fPt
1138 << "s2y=" << s2yTrk
1139 << "r0=" << rowId[0]
1140 << "dy0=" << &vdy[0]
1141 << "m0=" << mean[0]
1142 << "s0=" << syDis[0]
1143 << "r1=" << rowId[1]
1144 << "dy1=" << &vdy[1]
1145 << "m1=" << mean[1]
1146 << "s1=" << syDis[1]
1147 << "r2=" << rowId[2]
1148 << "dy2=" << &vdy[2]
1149 << "m2=" << mean[2]
1150 << "s2=" << syDis[2]
1151 << "\n";
1152 }
1153
1154
1155 // analyze gap in rows attached
1156 if(kRowSelection){
1157 SetErrorMsg(kAttachRowGap);
1158 Int_t rowRemove(-1);
1159 if(nr==2){ // select based on minimum distance to track projection
1160 if(TMath::Abs(mean[0])<TMath::Abs(mean[1])){
1161 if(nrow[1]>nrow[0]) AliDebug(2, Form("Conflicting mean[%f < %f] but ncl[%d < %d].", TMath::Abs(mean[0]), TMath::Abs(mean[1]), nrow[0], nrow[1]));
1162 }else{
1163 if(nrow[1]<nrow[0]) AliDebug(2, Form("Conflicting mean[%f > %f] but ncl[%d > %d].", TMath::Abs(mean[0]), TMath::Abs(mean[1]), nrow[0], nrow[1]));
1164 Swap(nrow[0],nrow[1]); Swap(rowId[0],rowId[1]);
1165 Swap(mean[0],mean[1]); Swap(syDis[0],syDis[1]);
1166 }
1167 rowRemove=1; nr=1;
1168 } else if(nr==3){ // select based on 2 consecutive rows
1169 if(rowId[1]==rowId[0]+1 && rowId[1]!=rowId[2]-1){
1170 nr=2;rowRemove=2;
1171 } else if(rowId[1]!=rowId[0]+1 && rowId[1]==rowId[2]-1){
1172 Swap(nrow[0],nrow[2]); Swap(rowId[0],rowId[2]);
1173 Swap(mean[0],mean[2]); Swap(syDis[0],syDis[2]);
1174 nr=2; rowRemove=2;
1175 }
29b87567 1176 }
560e5c05 1177 if(rowRemove>0){nrow[rowRemove]=0; rowId[rowRemove]=-1;}
29b87567 1178 }
560e5c05 1179 AliDebug(4, Form(" Ncl[%d[%d] + %d[%d] + %d[%d]]", nrow[0], rowId[0], nrow[1], rowId[1], nrow[2], rowId[2]));
1180
1181 if(nr==3){
1182 SetBit(kRowCross, kTRUE); // mark pad row crossing
7c3eecb8 1183 SetErrorMsg(kAttachRow);
560e5c05 1184 const Float_t am[]={TMath::Abs(mean[0]), TMath::Abs(mean[1]), TMath::Abs(mean[2])};
1185 AliDebug(4, Form("complex row configuration\n"
1186 " r[%d] n[%d] m[%6.3f] s[%6.3f]\n"
1187 " r[%d] n[%d] m[%6.3f] s[%6.3f]\n"
1188 " r[%d] n[%d] m[%6.3f] s[%6.3f]\n"
1189 , rowId[0], nrow[0], am[0], syDis[0]
1190 , rowId[1], nrow[1], am[1], syDis[1]
1191 , rowId[2], nrow[2], am[2], syDis[2]));
1192 Int_t id[]={0,1,2}; TMath::Sort(3, am, id, kFALSE);
1193 // backup
1194 Int_t rnn[3]; memcpy(rnn, nrow, 3*sizeof(Int_t));
1195 Int_t rid[3]; memcpy(rid, rowId, 3*sizeof(Int_t));
1196 Double_t rm[3]; memcpy(rm, mean, 3*sizeof(Double_t));
1197 Double_t rs[3]; memcpy(rs, syDis, 3*sizeof(Double_t));
1198 nrow[0]=rnn[id[0]]; rowId[0]=rid[id[0]]; mean[0]=rm[id[0]]; syDis[0]=rs[id[0]];
1199 nrow[1]=rnn[id[1]]; rowId[1]=rid[id[1]]; mean[1]=rm[id[1]]; syDis[1]=rs[id[1]];
1200 nrow[2]=0; rowId[2]=-1; mean[2] = 1.e3; syDis[2] = 1.e3;
1201 AliDebug(4, Form("solved configuration\n"
1202 " r[%d] n[%d] m[%+6.3f] s[%6.3f]\n"
1203 " r[%d] n[%d] m[%+6.3f] s[%6.3f]\n"
1204 " r[%d] n[%d] m[%+6.3f] s[%6.3f]\n"
1205 , rowId[0], nrow[0], mean[0], syDis[0]
1206 , rowId[1], nrow[1], mean[1], syDis[1]
1207 , rowId[2], nrow[2], mean[2], syDis[2]));
1208 nr=2;
1209 } else if(nr==2) {
1210 SetBit(kRowCross, kTRUE); // mark pad row crossing
1211 if(nrow[1] > nrow[0]){ // swap row order
1212 Swap(nrow[0],nrow[1]); Swap(rowId[0],rowId[1]);
1213 Swap(mean[0],mean[1]); Swap(syDis[0],syDis[1]);
1214 }
ee8fb199 1215 }
560e5c05 1216
b1957d3c 1217 // Select and store clusters
1218 // We should consider here :
1219 // 1. How far is the chamber boundary
1220 // 2. How big is the mean
560e5c05 1221 Int_t n(0); Float_t dyc[kNclusters]; memset(dyc,0,kNclusters*sizeof(Float_t));
b1957d3c 1222 for (Int_t ir = 0; ir < nr; ir++) {
560e5c05 1223 Int_t jr(rowId[ir]);
1224 AliDebug(4, Form(" Attaching Ncl[%d]=%d ...", jr, ncl[jr]));
b1957d3c 1225 for (Int_t ic = 0; ic < ncl[jr]; ic++) {
3044dfe5 1226 if(!blst[jr][ic])continue;
1227 c = clst[jr][ic];
560e5c05 1228 Int_t it(c->GetPadTime());
1229 Int_t idx(it+kNtb*ir);
6ad5e6b2 1230 if(fClusters[idx]){
560e5c05 1231 AliDebug(4, Form("Many cluster candidates on row[%2d] tb[%2d].", jr, it));
1232 // TODO should save also the information on where the multiplicity happened and its size
6ad5e6b2 1233 SetErrorMsg(kAttachMultipleCl);
560e5c05 1234 // TODO should also compare with mean and sigma for this row
1235 if(yres[jr][ic] > dyc[idx]) continue;
6ad5e6b2 1236 }
1237
b1957d3c 1238 // TODO proper indexing of clusters !!
6ad5e6b2 1239 fIndexes[idx] = chamber->GetTB(it)->GetGlobalIndex(idxs[jr][ic]);
1240 fClusters[idx] = c;
560e5c05 1241 dyc[idx] = yres[jr][ic];
3e778975 1242 n++;
b1957d3c 1243 }
560e5c05 1244 }
6ad5e6b2 1245 SetN(n);
b1957d3c 1246
29b87567 1247 // number of minimum numbers of clusters expected for the tracklet
6ad5e6b2 1248 if (GetN() < kClmin){
560e5c05 1249 AliDebug(1, Form("NOT ENOUGH CLUSTERS %d ATTACHED TO THE TRACKLET [min %d] FROM FOUND %d.", GetN(), kClmin, n));
7c3eecb8 1250 SetErrorMsg(kAttachClAttach);
e4f2f73d 1251 return kFALSE;
1252 }
0906e73e 1253
e3cf3d02 1254 // Load calibration parameters for this tracklet
1255 Calibrate();
b1957d3c 1256
1257 // calculate dx for time bins in the drift region (calibration aware)
a2abcbc5 1258 Float_t x[2] = {0.,0.}; Int_t tb[2]={0,0};
1259 for (Int_t it = t0, irp=0; irp<2 && it < AliTRDtrackerV1::GetNTimeBins(); it++) {
b1957d3c 1260 if(!fClusters[it]) continue;
1261 x[irp] = fClusters[it]->GetX();
a2abcbc5 1262 tb[irp] = fClusters[it]->GetLocalTimeBin();
b1957d3c 1263 irp++;
e3cf3d02 1264 }
d86ed84c 1265 Int_t dtb = tb[1] - tb[0];
1266 fdX = dtb ? (x[0] - x[1]) / dtb : 0.15;
29b87567 1267 return kTRUE;
e4f2f73d 1268}
1269
03cef9b2 1270//____________________________________________________________
1271void AliTRDseedV1::Bootstrap(const AliTRDReconstructor *rec)
1272{
1273// Fill in all derived information. It has to be called after recovery from file or HLT.
1274// The primitive data are
1275// - list of clusters
1276// - detector (as the detector will be removed from clusters)
1277// - position of anode wire (fX0) - temporary
1278// - track reference position and direction
1279// - momentum of the track
1280// - time bin length [cm]
1281//
1282// A.Bercuci <A.Bercuci@gsi.de> Oct 30th 2008
1283//
4d6aee34 1284 fkReconstructor = rec;
03cef9b2 1285 AliTRDgeometry g;
2eb10c34 1286 SetPadPlane(g.GetPadPlane(fDet));
1287
e3cf3d02 1288 //fSnp = fYref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);
1289 //fTgl = fZref[1];
3e778975 1290 Int_t n = 0, nshare = 0, nused = 0;
03cef9b2 1291 AliTRDcluster **cit = &fClusters[0];
8d2bec9e 1292 for(Int_t ic = kNclusters; ic--; cit++){
03cef9b2 1293 if(!(*cit)) return;
3e778975 1294 n++;
1295 if((*cit)->IsShared()) nshare++;
1296 if((*cit)->IsUsed()) nused++;
03cef9b2 1297 }
3e778975 1298 SetN(n); SetNUsed(nused); SetNShared(nshare);
e3cf3d02 1299 Fit();
03cef9b2 1300 CookLabels();
1301 GetProbability();
1302}
1303
1304
e4f2f73d 1305//____________________________________________________________________
2eb10c34 1306Bool_t AliTRDseedV1::Fit(UChar_t opt)
e4f2f73d 1307{
16cca13f 1308//
1309// Linear fit of the clusters attached to the tracklet
1310//
1311// Parameters :
2eb10c34 1312// - opt : switch for tilt pad correction of cluster y position. Options are
1313// 0 no correction [default]
1314// 1 full tilt correction [dz/dx and z0]
1315// 2 pseudo tilt correction [dz/dx from pad-chamber geometry]
1316//
16cca13f 1317// Output :
1318// True if successful
1319//
1320// Detailed description
1321//
1322// Fit in the xy plane
1323//
1fd9389f 1324// The fit is performed to estimate the y position of the tracklet and the track
1325// angle in the bending plane. The clusters are represented in the chamber coordinate
1326// system (with respect to the anode wire - see AliTRDtrackerV1::FollowBackProlongation()
1327// on how this is set). The x and y position of the cluster and also their variances
1328// are known from clusterizer level (see AliTRDcluster::GetXloc(), AliTRDcluster::GetYloc(),
1329// AliTRDcluster::GetSX() and AliTRDcluster::GetSY()).
1330// If gaussian approximation is used to calculate y coordinate of the cluster the position
1331// is recalculated taking into account the track angle. The general formula to calculate the
1332// error of cluster position in the gaussian approximation taking into account diffusion and track
1333// inclination is given for TRD by:
1334// BEGIN_LATEX
1335// #sigma^{2}_{y} = #sigma^{2}_{PRF} + #frac{x#delta_{t}^{2}}{(1+tg(#alpha_{L}))^{2}} + #frac{x^{2}tg^{2}(#phi-#alpha_{L})tg^{2}(#alpha_{L})}{12}
1336// END_LATEX
1337//
1338// Since errors are calculated only in the y directions, radial errors (x direction) are mapped to y
1339// by projection i.e.
1340// BEGIN_LATEX
1341// #sigma_{x|y} = tg(#phi) #sigma_{x}
1342// END_LATEX
1343// and also by the lorentz angle correction
1344//
1345// Fit in the xz plane
1346//
1347// The "fit" is performed to estimate the radial position (x direction) where pad row cross happens.
1348// If no pad row crossing the z position is taken from geometry and radial position is taken from the xy
1349// fit (see below).
1350//
1351// There are two methods to estimate the radial position of the pad row cross:
1352// 1. leading cluster radial position : Here the lower part of the tracklet is considered and the last
1353// cluster registered (at radial x0) on this segment is chosen to mark the pad row crossing. The error
1354// of the z estimate is given by :
1355// BEGIN_LATEX
1356// #sigma_{z} = tg(#theta) #Delta x_{x_{0}}/12
1357// END_LATEX
1358// The systematic errors for this estimation are generated by the following sources:
1359// - no charge sharing between pad rows is considered (sharp cross)
1360// - missing cluster at row cross (noise peak-up, under-threshold signal etc.).
1361//
1362// 2. charge fit over the crossing point : Here the full energy deposit along the tracklet is considered
1363// to estimate the position of the crossing by a fit in the qx plane. The errors in the q directions are
1364// parameterized as s_q = q^2. The systematic errors for this estimation are generated by the following sources:
1365// - no general model for the qx dependence
1366// - physical fluctuations of the charge deposit
1367// - gain calibration dependence
1368//
1369// Estimation of the radial position of the tracklet
16cca13f 1370//
1fd9389f 1371// For pad row cross the radial position is taken from the xz fit (see above). Otherwise it is taken as the
1372// interpolation point of the tracklet i.e. the point where the error in y of the fit is minimum. The error
1373// in the y direction of the tracklet is (see AliTRDseedV1::GetCovAt()):
1374// BEGIN_LATEX
1375// #sigma_{y} = #sigma^{2}_{y_{0}} + 2xcov(y_{0}, dy/dx) + #sigma^{2}_{dy/dx}
1376// END_LATEX
1377// and thus the radial position is:
1378// BEGIN_LATEX
1379// x = - cov(y_{0}, dy/dx)/#sigma^{2}_{dy/dx}
1380// END_LATEX
1381//
1382// Estimation of tracklet position error
1383//
1384// The error in y direction is the error of the linear fit at the radial position of the tracklet while in the z
1385// direction is given by the cluster error or pad row cross error. In case of no pad row cross this is given by:
1386// BEGIN_LATEX
1387// #sigma_{y} = #sigma^{2}_{y_{0}} - 2cov^{2}(y_{0}, dy/dx)/#sigma^{2}_{dy/dx} + #sigma^{2}_{dy/dx}
1388// #sigma_{z} = Pad_{length}/12
1389// END_LATEX
1390// For pad row cross the full error is calculated at the radial position of the crossing (see above) and the error
1391// in z by the width of the crossing region - being a matter of parameterization.
1392// BEGIN_LATEX
1393// #sigma_{z} = tg(#theta) #Delta x_{x_{0}}/12
1394// END_LATEX
1395// In case of no tilt correction (default in the barrel tracking) the tilt is taken into account by the rotation of
1396// the covariance matrix. See AliTRDseedV1::GetCovAt() for details.
1397//
1398// Author
1399// A.Bercuci <A.Bercuci@gsi.de>
e4f2f73d 1400
a723055f 1401 if(!fkReconstructor){
1402 AliError("The tracklet needs the reconstruction setup. Please initialize by SetReconstructor().");
1403 return kFALSE;
1404 }
b72f4eaf 1405 if(!IsCalibrated()) Calibrate();
2eb10c34 1406 if(opt>2){
7e5954f0 1407 AliWarning(Form("Option [%d] outside range [0, 2]. Using default",opt));
2eb10c34 1408 opt=0;
1409 }
e3cf3d02 1410
29b87567 1411 const Int_t kClmin = 8;
2eb10c34 1412 const Float_t kScalePulls = 10.; // factor to scale y pulls - NOT UNDERSTOOD
2f7d6ac8 1413 // get track direction
1414 Double_t y0 = fYref[0];
1415 Double_t dydx = fYref[1];
1416 Double_t z0 = fZref[0];
1417 Double_t dzdx = fZref[1];
ae4e8b84 1418
5f1ae1e7 1419 AliTRDtrackerV1::AliTRDLeastSquare fitterY;
1420 AliTRDtrackerV1::AliTRDLeastSquare fitterZ;
f301a656 1421
29b87567 1422 // book cluster information
8d2bec9e 1423 Double_t qc[kNclusters], xc[kNclusters], yc[kNclusters], zc[kNclusters], sy[kNclusters];
e3cf3d02 1424
2eb10c34 1425 Bool_t tilt(opt==1) // full tilt correction
1426 ,pseudo(opt==2) // pseudo tilt correction
1427 ,rc(IsRowCross()) // row cross candidate
1428 ,kDZDX(IsPrimary());// switch dzdx calculation for barrel primary tracks
1429 Int_t n(0); // clusters used in fit
1430 AliTRDcluster *c(NULL), *cc(NULL), **jc = &fClusters[0];
fc0882f3 1431 const AliTRDrecoParam* const recoParam = fkReconstructor->GetRecoParam(); //the dynamic cast in GetRecoParam is slow, so caching the pointer to it
2eb10c34 1432
1433 const Char_t *tcName[]={"NONE", "FULL", "HALF"};
1434 AliDebug(2, Form("Options : TC[%s] dzdx[%c]", tcName[opt], kDZDX?'Y':'N'));
1435
1436 for (Int_t ic=0; ic<kNclusters; ic++, ++jc) {
1437 xc[ic] = -1.; yc[ic] = 999.; zc[ic] = 999.; sy[ic] = 0.;
9eb2d46c 1438 if(!(c = (*jc))) continue;
29b87567 1439 if(!c->IsInChamber()) continue;
2eb10c34 1440 // compute pseudo tilt correction
1441 if(kDZDX){
1442 fZfit[0] = c->GetZ();
1443 if(rc){
1444 for(Int_t kc=AliTRDseedV1::kNtb; kc<AliTRDseedV1::kNclusters; kc++){
1445 if(!(cc=fClusters[kc])) continue;
1446 if(!cc->IsInChamber()) continue;
1447 fZfit[0] += cc->GetZ(); fZfit[0] *= 0.5;
1448 break;
1449 }
1450 }
1451 fZfit[1] = fZfit[0]/fX0;
1452 if(rc){
1453 fZfit[0] += fZfit[1]*0.5*AliTRDgeometry::CdrHght();
1454 fZfit[1] = fZfit[0]/fX0;
1455 }
1456 kDZDX=kFALSE;
1457 }
9462866a 1458
29b87567 1459 Float_t w = 1.;
1460 if(c->GetNPads()>4) w = .5;
1461 if(c->GetNPads()>5) w = .2;
010d62b0 1462
1fd9389f 1463 // cluster charge
dd8059a8 1464 qc[n] = TMath::Abs(c->GetQ());
1fd9389f 1465 // pad row of leading
1466
b72f4eaf 1467 xc[n] = fX0 - c->GetX();
1468
1fd9389f 1469 // Recalculate cluster error based on tracking information
2eb10c34 1470 c->SetSigmaY2(fS2PRF, fDiffT, fExB, xc[n], -1./*zcorr?zt:-1.*/, dydx);
c79857d5 1471 c->SetSigmaZ2(fPad[0]*fPad[0]/12.); // for HLT
1fd9389f 1472 sy[n] = TMath::Sqrt(c->GetSigmaY2());
1473
fc0882f3 1474 yc[n] = recoParam->UseGAUS() ?
1fd9389f 1475 c->GetYloc(y0, sy[n], GetPadWidth()): c->GetY();
1476 zc[n] = c->GetZ();
2eb10c34 1477
1478 //optional r-phi correction
1479 //printf(" n[%2d] yc[%7.5f] ", n, yc[n]);
1480 Float_t correction(0.);
1481 if(tilt) correction = fPad[2]*(xc[n]*dzdx + zc[n] - z0);
1482 else if(pseudo) correction = fPad[2]*(xc[n]*fZfit[1] + zc[n]-fZfit[0]);
1483 yc[n]-=correction;
1484 //printf("corr(%s%s)[%7.5f] yc1[%7.5f]\n", (tilt?"TC":""), (zcorr?"PC":""), correction, yc[n]);
1fd9389f 1485
fbe11be7 1486 AliDebug(5, Form(" tb[%2d] dx[%6.3f] y[%6.2f+-%6.3f]", c->GetLocalTimeBin(), xc[n], yc[n], sy[n]));
903326c1 1487 fitterY.AddPoint(&xc[n], yc[n], sy[n]);
2eb10c34 1488 if(rc) fitterZ.AddPoint(&xc[n], qc[n]*(ic<kNtb?1.:-1.), 1.);
dd8059a8 1489 n++;
29b87567 1490 }
3044dfe5 1491
47d5d320 1492 // to few clusters
c79857d5 1493 if (n < kClmin){
c388cdcb 1494 AliDebug(1, Form("Not enough clusters to fit. Clusters: Attached[%d] Fit[%d].", GetN(), n));
2eb10c34 1495 SetErrorMsg(kFitCl);
c79857d5 1496 return kFALSE;
1497 }
d937ad7a 1498 // fit XY
903326c1 1499 if(!fitterY.Eval()){
c388cdcb 1500 AliDebug(1, "Fit Y failed.");
2eb10c34 1501 SetErrorMsg(kFitFailedY);
903326c1 1502 return kFALSE;
1503 }
5f1ae1e7 1504 fYfit[0] = fitterY.GetFunctionParameter(0);
1505 fYfit[1] = -fitterY.GetFunctionParameter(1);
d937ad7a 1506 // store covariance
5f1ae1e7 1507 Double_t p[3];
1508 fitterY.GetCovarianceMatrix(p);
2eb10c34 1509 fCov[0] = kScalePulls*p[1]; // variance of y0
1510 fCov[1] = kScalePulls*p[2]; // covariance of y0, dydx
1511 fCov[2] = kScalePulls*p[0]; // variance of dydx
b1957d3c 1512 // the ref radial position is set at the minimum of
1513 // the y variance of the tracklet
b72f4eaf 1514 fX = -fCov[1]/fCov[2];
2eb10c34 1515 fS2Y = fCov[0] +2.*fX*fCov[1] + fX*fX*fCov[2];
1516
903326c1 1517 Float_t xs=fX+.5*AliTRDgeometry::CamHght();
1518 if(xs < 0. || xs > AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght()){
1519 AliDebug(1, Form("Ref radial position ouside chamber x[%5.2f].", fX));
2eb10c34 1520 SetErrorMsg(kFitFailedY);
903326c1 1521 return kFALSE;
1522 }
b1957d3c 1523
2eb10c34 1524/* // THE LEADING CLUSTER METHOD for z fit
1fd9389f 1525 Float_t xMin = fX0;
b72f4eaf 1526 Int_t ic=n=kNclusters-1; jc = &fClusters[ic];
1fd9389f 1527 AliTRDcluster *c0 =0x0, **kc = &fClusters[kNtb-1];
1528 for(; ic>kNtb; ic--, --jc, --kc){
1529 if((c0 = (*kc)) && c0->IsInChamber() && (xMin>c0->GetX())) xMin = c0->GetX();
1530 if(!(c = (*jc))) continue;
1531 if(!c->IsInChamber()) continue;
1532 zc[kNclusters-1] = c->GetZ();
1533 fX = fX0 - c->GetX();
1534 }
1535 fZfit[0] = .5*(zc[0]+zc[kNclusters-1]); fZfit[1] = 0.;
1536 // Error parameterization
1537 fS2Z = fdX*fZref[1];
e355f67a 1538 fS2Z *= fS2Z; fS2Z *= 0.2887; // 1/sqrt(12)*/
1539
2eb10c34 1540 // fit QZ
1541 if(opt!=1 && IsRowCross()){
1542 if(!fitterZ.Eval()) SetErrorMsg(kFitFailedZ);
4ecadb52 1543 if(!HasError(kFitFailedZ) && TMath::Abs(fitterZ.GetFunctionParameter(1))>1.e-10){
2eb10c34 1544 // TODO - one has to recalculate xy fit based on
1545 // better knowledge of z position
1546// Double_t x = -fitterZ.GetFunctionParameter(0)/fitterZ.GetFunctionParameter(1);
1547// Double_t z0 = .5*(zc[0]+zc[n-1]);
1548// fZfit[0] = z0 + fZfit[1]*x;
1549// fZfit[1] = fZfit[0]/fX0;
1550// redo fit on xy plane
b72f4eaf 1551 }
c850c351 1552 // temporary external error parameterization
1553 fS2Z = 0.05+0.4*TMath::Abs(fZref[1]); fS2Z *= fS2Z;
1554 // TODO correct formula
1555 //fS2Z = sigma_x*TMath::Abs(fZref[1]);
b1957d3c 1556 } else {
2eb10c34 1557 //fZfit[0] = zc[0] + dzdx*0.5*AliTRDgeometry::CdrHght();
dd8059a8 1558 fS2Z = GetPadLength()*GetPadLength()/12.;
29b87567 1559 }
29b87567 1560 return kTRUE;
e4f2f73d 1561}
1562
e4f2f73d 1563
f29f13a6 1564/*
e3cf3d02 1565//_____________________________________________________________________________
1566void AliTRDseedV1::FitMI()
1567{
1568//
1569// Fit the seed.
1570// Marian Ivanov's version
1571//
1572// linear fit on the y direction with respect to the reference direction.
1573// The residuals for each x (x = xc - x0) are deduced from:
1574// dy = y - yt (1)
1575// the tilting correction is written :
1576// y = yc + h*(zc-zt) (2)
1577// yt = y0+dy/dx*x (3)
1578// zt = z0+dz/dx*x (4)
1579// from (1),(2),(3) and (4)
1580// dy = yc - y0 - (dy/dx + h*dz/dx)*x + h*(zc-z0)
1581// the last term introduces the correction on y direction due to tilting pads. There are 2 ways to account for this:
1582// 1. use tilting correction for calculating the y
1583// 2. neglect tilting correction here and account for it in the error parametrization of the tracklet.
1584 const Float_t kRatio = 0.8;
1585 const Int_t kClmin = 5;
1586 const Float_t kmaxtan = 2;
1587
1588 if (TMath::Abs(fYref[1]) > kmaxtan){
1589 //printf("Exit: Abs(fYref[1]) = %3.3f, kmaxtan = %3.3f\n", TMath::Abs(fYref[1]), kmaxtan);
1590 return; // Track inclined too much
1591 }
1592
1593 Float_t sigmaexp = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction
dd8059a8 1594 Float_t ycrosscor = GetPadLength() * GetTilt() * 0.5; // Y correction for crossing
e3cf3d02 1595 Int_t fNChange = 0;
1596
1597 Double_t sumw;
1598 Double_t sumwx;
1599 Double_t sumwx2;
1600 Double_t sumwy;
1601 Double_t sumwxy;
1602 Double_t sumwz;
1603 Double_t sumwxz;
1604
1605 // Buffering: Leave it constant fot Performance issues
1606 Int_t zints[kNtb]; // Histograming of the z coordinate
1607 // Get 1 and second max probable coodinates in z
1608 Int_t zouts[2*kNtb];
1609 Float_t allowedz[kNtb]; // Allowed z for given time bin
1610 Float_t yres[kNtb]; // Residuals from reference
dd8059a8 1611 //Float_t anglecor = GetTilt() * fZref[1]; // Correction to the angle
e3cf3d02 1612
1613 Float_t pos[3*kNtb]; memset(pos, 0, 3*kNtb*sizeof(Float_t));
1614 Float_t *fX = &pos[0], *fY = &pos[kNtb], *fZ = &pos[2*kNtb];
1615
1616 Int_t fN = 0; AliTRDcluster *c = 0x0;
1617 fN2 = 0;
1618 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1619 yres[i] = 10000.0;
1620 if (!(c = fClusters[i])) continue;
1621 if(!c->IsInChamber()) continue;
1622 // Residual y
dd8059a8 1623 //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] + GetTilt()*(fZ[i] - fZref[0]);
e3cf3d02 1624 fX[i] = fX0 - c->GetX();
1625 fY[i] = c->GetY();
1626 fZ[i] = c->GetZ();
dd8059a8 1627 yres[i] = fY[i] - GetTilt()*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
e3cf3d02 1628 zints[fN] = Int_t(fZ[i]);
1629 fN++;
1630 }
1631
1632 if (fN < kClmin){
1633 //printf("Exit fN < kClmin: fN = %d\n", fN);
1634 return;
1635 }
1636 Int_t nz = AliTRDtrackerV1::Freq(fN, zints, zouts, kFALSE);
1637 Float_t fZProb = zouts[0];
1638 if (nz <= 1) zouts[3] = 0;
1639 if (zouts[1] + zouts[3] < kClmin) {
1640 //printf("Exit zouts[1] = %d, zouts[3] = %d\n",zouts[1],zouts[3]);
1641 return;
1642 }
1643
1644 // Z distance bigger than pad - length
1645 if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) zouts[3] = 0;
1646
1647 Int_t breaktime = -1;
1648 Bool_t mbefore = kFALSE;
1649 Int_t cumul[kNtb][2];
1650 Int_t counts[2] = { 0, 0 };
1651
1652 if (zouts[3] >= 3) {
1653
1654 //
1655 // Find the break time allowing one chage on pad-rows
1656 // with maximal number of accepted clusters
1657 //
1658 fNChange = 1;
1659 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1660 cumul[i][0] = counts[0];
1661 cumul[i][1] = counts[1];
1662 if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++;
1663 if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++;
1664 }
1665 Int_t maxcount = 0;
1666 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1667 Int_t after = cumul[AliTRDtrackerV1::GetNTimeBins()][0] - cumul[i][0];
1668 Int_t before = cumul[i][1];
1669 if (after + before > maxcount) {
1670 maxcount = after + before;
1671 breaktime = i;
1672 mbefore = kFALSE;
1673 }
1674 after = cumul[AliTRDtrackerV1::GetNTimeBins()-1][1] - cumul[i][1];
1675 before = cumul[i][0];
1676 if (after + before > maxcount) {
1677 maxcount = after + before;
1678 breaktime = i;
1679 mbefore = kTRUE;
1680 }
1681 }
1682 breaktime -= 1;
1683 }
1684
1685 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1686 if (i > breaktime) allowedz[i] = mbefore ? zouts[2] : zouts[0];
1687 if (i <= breaktime) allowedz[i] = (!mbefore) ? zouts[2] : zouts[0];
1688 }
1689
1690 if (((allowedz[0] > allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] < 0)) ||
1691 ((allowedz[0] < allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] > 0))) {
1692 //
1693 // Tracklet z-direction not in correspondance with track z direction
1694 //
1695 fNChange = 0;
1696 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1697 allowedz[i] = zouts[0]; // Only longest taken
1698 }
1699 }
1700
1701 if (fNChange > 0) {
1702 //
1703 // Cross pad -row tracklet - take the step change into account
1704 //
1705 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1706 if (!fClusters[i]) continue;
1707 if(!fClusters[i]->IsInChamber()) continue;
1708 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
1709 // Residual y
dd8059a8 1710 //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] + GetTilt()*(fZ[i] - fZref[0]);
1711 yres[i] = fY[i] - GetTilt()*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
f29f13a6 1712// if (TMath::Abs(fZ[i] - fZProb) > 2) {
dd8059a8 1713// if (fZ[i] > fZProb) yres[i] += GetTilt() * GetPadLength();
1714// if (fZ[i] < fZProb) yres[i] -= GetTilt() * GetPadLength();
f29f13a6 1715 }
e3cf3d02 1716 }
1717 }
1718
1719 Double_t yres2[kNtb];
1720 Double_t mean;
1721 Double_t sigma;
1722 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1723 if (!fClusters[i]) continue;
1724 if(!fClusters[i]->IsInChamber()) continue;
1725 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
1726 yres2[fN2] = yres[i];
1727 fN2++;
1728 }
1729 if (fN2 < kClmin) {
1730 //printf("Exit fN2 < kClmin: fN2 = %d\n", fN2);
1731 fN2 = 0;
1732 return;
1733 }
1734 AliMathBase::EvaluateUni(fN2,yres2,mean,sigma, Int_t(fN2*kRatio-2.));
1735 if (sigma < sigmaexp * 0.8) {
1736 sigma = sigmaexp;
1737 }
1738 //Float_t fSigmaY = sigma;
1739
1740 // Reset sums
1741 sumw = 0;
1742 sumwx = 0;
1743 sumwx2 = 0;
1744 sumwy = 0;
1745 sumwxy = 0;
1746 sumwz = 0;
1747 sumwxz = 0;
1748
1749 fN2 = 0;
1750 Float_t fMeanz = 0;
1751 Float_t fMPads = 0;
1752 fUsable = 0;
1753 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1754 if (!fClusters[i]) continue;
1755 if (!fClusters[i]->IsInChamber()) continue;
1756 if (TMath::Abs(fZ[i] - allowedz[i]) > 2){fClusters[i] = 0x0; continue;}
1757 if (TMath::Abs(yres[i] - mean) > 4.0 * sigma){fClusters[i] = 0x0; continue;}
1758 SETBIT(fUsable,i);
1759 fN2++;
1760 fMPads += fClusters[i]->GetNPads();
1761 Float_t weight = 1.0;
1762 if (fClusters[i]->GetNPads() > 4) weight = 0.5;
1763 if (fClusters[i]->GetNPads() > 5) weight = 0.2;
1764
1765
1766 Double_t x = fX[i];
1767 //printf("x = %7.3f dy = %7.3f fit %7.3f\n", x, yres[i], fY[i]-yres[i]);
1768
1769 sumw += weight;
1770 sumwx += x * weight;
1771 sumwx2 += x*x * weight;
1772 sumwy += weight * yres[i];
1773 sumwxy += weight * (yres[i]) * x;
1774 sumwz += weight * fZ[i];
1775 sumwxz += weight * fZ[i] * x;
1776
1777 }
1778
1779 if (fN2 < kClmin){
1780 //printf("Exit fN2 < kClmin(2): fN2 = %d\n",fN2);
1781 fN2 = 0;
1782 return;
1783 }
1784 fMeanz = sumwz / sumw;
1785 Float_t correction = 0;
1786 if (fNChange > 0) {
1787 // Tracklet on boundary
1788 if (fMeanz < fZProb) correction = ycrosscor;
1789 if (fMeanz > fZProb) correction = -ycrosscor;
1790 }
1791
1792 Double_t det = sumw * sumwx2 - sumwx * sumwx;
1793 fYfit[0] = (sumwx2 * sumwy - sumwx * sumwxy) / det;
1794 fYfit[1] = (sumw * sumwxy - sumwx * sumwy) / det;
1795
1796 fS2Y = 0;
1797 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1798 if (!TESTBIT(fUsable,i)) continue;
1799 Float_t delta = yres[i] - fYfit[0] - fYfit[1] * fX[i];
1800 fS2Y += delta*delta;
1801 }
1802 fS2Y = TMath::Sqrt(fS2Y / Float_t(fN2-2));
1803 // TEMPORARY UNTIL covariance properly calculated
1804 fS2Y = TMath::Max(fS2Y, Float_t(.1));
1805
1806 fZfit[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det;
1807 fZfit[1] = (sumw * sumwxz - sumwx * sumwz) / det;
1808// fYfitR[0] += fYref[0] + correction;
1809// fYfitR[1] += fYref[1];
1810// fYfit[0] = fYfitR[0];
1811 fYfit[1] = -fYfit[1];
1812
1813 UpdateUsed();
f29f13a6 1814}*/
e3cf3d02 1815
e4f2f73d 1816//___________________________________________________________________
203967fc 1817void AliTRDseedV1::Print(Option_t *o) const
e4f2f73d 1818{
1819 //
1820 // Printing the seedstatus
1821 //
1822
b72f4eaf 1823 AliInfo(Form("Det[%3d] X0[%7.2f] Pad{L[%5.2f] W[%5.2f] Tilt[%+6.2f]}", fDet, fX0, GetPadLength(), GetPadWidth(), GetTilt()));
dd8059a8 1824 AliInfo(Form("N[%2d] Nused[%2d] Nshared[%2d] [%d]", GetN(), GetNUsed(), GetNShared(), fN));
b72f4eaf 1825 AliInfo(Form("FLAGS : RC[%c] Kink[%c] SA[%c]", IsRowCross()?'y':'n', IsKink()?'y':'n', IsStandAlone()?'y':'n'));
525f399b 1826 AliInfo(Form("CALIB PARAMS : T0[%5.2f] Vd[%5.2f] s2PRF[%5.2f] ExB[%5.2f] Dl[%5.2f] Dt[%5.2f]", fT0, fVD, fS2PRF, fExB, fDiffL, fDiffT));
dd8059a8 1827
1828 Double_t cov[3], x=GetX();
1829 GetCovAt(x, cov);
1830 AliInfo(" | x[cm] | y[cm] | z[cm] | dydx | dzdx |");
1831 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]));
16cca13f 1832 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[0]), fZref[0]-fX*fYref[1], TMath::Sqrt(fRefCov[2]), fYref[1], fZref[1]))
ee8fb199 1833 AliInfo(Form("P / Pt [GeV/c] = %f / %f", GetMomentum(), fPt));
68f9b6bd 1834 if(IsStandAlone()) AliInfo(Form("C Rieman / Vertex [1/cm] = %f / %f", fC[0], fC[1]));
ee8fb199 1835 AliInfo(Form("dEdx [a.u.] = %f / %f / %f / %f / %f/ %f / %f / %f", fdEdx[0], fdEdx[1], fdEdx[2], fdEdx[3], fdEdx[4], fdEdx[5], fdEdx[6], fdEdx[7]));
1836 AliInfo(Form("PID = %5.3f / %5.3f / %5.3f / %5.3f / %5.3f", fProb[0], fProb[1], fProb[2], fProb[3], fProb[4]));
203967fc 1837
1838 if(strcmp(o, "a")!=0) return;
1839
4dc4dc2e 1840 AliTRDcluster* const* jc = &fClusters[0];
8d2bec9e 1841 for(int ic=0; ic<kNclusters; ic++, jc++) {
4dc4dc2e 1842 if(!(*jc)) continue;
203967fc 1843 (*jc)->Print(o);
4dc4dc2e 1844 }
e4f2f73d 1845}
47d5d320 1846
203967fc 1847
1848//___________________________________________________________________
1849Bool_t AliTRDseedV1::IsEqual(const TObject *o) const
1850{
1851 // Checks if current instance of the class has the same essential members
1852 // as the given one
1853
1854 if(!o) return kFALSE;
1855 const AliTRDseedV1 *inTracklet = dynamic_cast<const AliTRDseedV1*>(o);
1856 if(!inTracklet) return kFALSE;
1857
1858 for (Int_t i = 0; i < 2; i++){
e3cf3d02 1859 if ( fYref[i] != inTracklet->fYref[i] ) return kFALSE;
1860 if ( fZref[i] != inTracklet->fZref[i] ) return kFALSE;
203967fc 1861 }
1862
4ecadb52 1863 if ( TMath::Abs(fS2Y - inTracklet->fS2Y)>1.e-10 ) return kFALSE;
1864 if ( TMath::Abs(GetTilt() - inTracklet->GetTilt())>1.e-10 ) return kFALSE;
1865 if ( TMath::Abs(GetPadLength() - inTracklet->GetPadLength())>1.e-10 ) return kFALSE;
203967fc 1866
8d2bec9e 1867 for (Int_t i = 0; i < kNclusters; i++){
e3cf3d02 1868// if ( fX[i] != inTracklet->GetX(i) ) return kFALSE;
1869// if ( fY[i] != inTracklet->GetY(i) ) return kFALSE;
1870// if ( fZ[i] != inTracklet->GetZ(i) ) return kFALSE;
1871 if ( fIndexes[i] != inTracklet->fIndexes[i] ) return kFALSE;
203967fc 1872 }
f29f13a6 1873// if ( fUsable != inTracklet->fUsable ) return kFALSE;
203967fc 1874
1875 for (Int_t i=0; i < 2; i++){
e3cf3d02 1876 if ( fYfit[i] != inTracklet->fYfit[i] ) return kFALSE;
1877 if ( fZfit[i] != inTracklet->fZfit[i] ) return kFALSE;
1878 if ( fLabels[i] != inTracklet->fLabels[i] ) return kFALSE;
203967fc 1879 }
1880
e3cf3d02 1881/* if ( fMeanz != inTracklet->GetMeanz() ) return kFALSE;
1882 if ( fZProb != inTracklet->GetZProb() ) return kFALSE;*/
3e778975 1883 if ( fN != inTracklet->fN ) return kFALSE;
1884 //if ( fNUsed != inTracklet->fNUsed ) return kFALSE;
e3cf3d02 1885 //if ( fFreq != inTracklet->GetFreq() ) return kFALSE;
1886 //if ( fNChange != inTracklet->GetNChange() ) return kFALSE;
203967fc 1887
4ecadb52 1888 if ( TMath::Abs(fC[0] - inTracklet->fC[0])>1.e-10 ) return kFALSE;
e3cf3d02 1889 //if ( fCC != inTracklet->GetCC() ) return kFALSE;
4ecadb52 1890 if ( TMath::Abs(fChi2 - inTracklet->fChi2)>1.e-10 ) return kFALSE;
203967fc 1891 // if ( fChi2Z != inTracklet->GetChi2Z() ) return kFALSE;
1892
e3cf3d02 1893 if ( fDet != inTracklet->fDet ) return kFALSE;
4ecadb52 1894 if ( TMath::Abs(fPt - inTracklet->fPt)>1.e-10 ) return kFALSE;
1895 if ( TMath::Abs(fdX - inTracklet->fdX)>1.e-10 ) return kFALSE;
203967fc 1896
8d2bec9e 1897 for (Int_t iCluster = 0; iCluster < kNclusters; iCluster++){
203967fc 1898 AliTRDcluster *curCluster = fClusters[iCluster];
e3cf3d02 1899 AliTRDcluster *inCluster = inTracklet->fClusters[iCluster];
203967fc 1900 if (curCluster && inCluster){
1901 if (! curCluster->IsEqual(inCluster) ) {
1902 curCluster->Print();
1903 inCluster->Print();
1904 return kFALSE;
1905 }
1906 } else {
1907 // if one cluster exists, and corresponding
1908 // in other tracklet doesn't - return kFALSE
1909 if(curCluster || inCluster) return kFALSE;
1910 }
1911 }
1912 return kTRUE;
1913}
5d401b45 1914