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