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