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