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