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