xxx
[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"
9dcc64cc 39#include "TTreeStream.h"
40#include "TGraphErrors.h"
e4f2f73d 41
42#include "AliLog.h"
43#include "AliMathBase.h"
4ecadb52 44#include "AliRieman.h"
d937ad7a 45#include "AliCDBManager.h"
e4f2f73d 46
4ecadb52 47#include "AliTRDReconstructor.h"
03cef9b2 48#include "AliTRDpadPlane.h"
9dcc64cc 49#include "AliTRDtransform.h"
e4f2f73d 50#include "AliTRDcluster.h"
f3d3af1b 51#include "AliTRDseedV1.h"
52#include "AliTRDtrackV1.h"
e4f2f73d 53#include "AliTRDcalibDB.h"
eb38ed55 54#include "AliTRDchamberTimeBin.h"
55#include "AliTRDtrackingChamber.h"
56#include "AliTRDtrackerV1.h"
e4f2f73d 57#include "AliTRDrecoParam.h"
a076fc2f 58#include "AliTRDCommonParam.h"
9dcc64cc 59#include "AliTRDtrackletOflHelper.h"
d937ad7a 60
9dcc64cc 61#include "Cal/AliTRDCalTrkAttach.h"
0906e73e 62#include "Cal/AliTRDCalPID.h"
d937ad7a 63#include "Cal/AliTRDCalROC.h"
64#include "Cal/AliTRDCalDet.h"
e4f2f73d 65
4ecadb52 66class AliTracker;
67
e4f2f73d 68ClassImp(AliTRDseedV1)
69
70//____________________________________________________________________
ae4e8b84 71AliTRDseedV1::AliTRDseedV1(Int_t det)
3e778975 72 :AliTRDtrackletBase()
4d6aee34 73 ,fkReconstructor(NULL)
74 ,fClusterIter(NULL)
e3cf3d02 75 ,fExB(0.)
76 ,fVD(0.)
77 ,fT0(0.)
78 ,fS2PRF(0.)
79 ,fDiffL(0.)
80 ,fDiffT(0.)
ae4e8b84 81 ,fClusterIdx(0)
7c3eecb8 82 ,fErrorMsg(0)
3e778975 83 ,fN(0)
ae4e8b84 84 ,fDet(det)
b25a5e9e 85 ,fPt(0.)
bcb6fb78 86 ,fdX(0.)
e3cf3d02 87 ,fX0(0.)
88 ,fX(0.)
89 ,fY(0.)
90 ,fZ(0.)
91 ,fS2Y(0.)
92 ,fS2Z(0.)
e3cf3d02 93 ,fChi2(0.)
e4f2f73d 94{
95 //
96 // Constructor
97 //
f301a656 98 memset(fIndexes,0xFF,kNclusters*sizeof(fIndexes[0]));
8d2bec9e 99 memset(fClusters, 0, kNclusters*sizeof(AliTRDcluster*));
2eb10c34 100 memset(fPad, 0, 4*sizeof(Float_t));
e3cf3d02 101 fYref[0] = 0.; fYref[1] = 0.;
102 fZref[0] = 0.; fZref[1] = 0.;
103 fYfit[0] = 0.; fYfit[1] = 0.;
104 fZfit[0] = 0.; fZfit[1] = 0.;
1f97f376 105 memset(fdEdx, 0, kNdEdxSlices*sizeof(Float_t));
29b87567 106 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = -1.;
e3cf3d02 107 fLabels[0]=-1; fLabels[1]=-1; // most freq MC labels
108 fLabels[2]=0; // number of different labels for tracklet
16cca13f 109 memset(fRefCov, 0, 7*sizeof(Double_t));
68f9b6bd 110 // stand alone curvature
111 fC[0] = 0.; fC[1] = 0.;
d937ad7a 112 // covariance matrix [diagonal]
113 // default sy = 200um and sz = 2.3 cm
114 fCov[0] = 4.e-4; fCov[1] = 0.; fCov[2] = 5.3;
f29f13a6 115 SetStandAlone(kFALSE);
e4f2f73d 116}
117
118//____________________________________________________________________
0906e73e 119AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
3e778975 120 :AliTRDtrackletBase((AliTRDtrackletBase&)ref)
4d6aee34 121 ,fkReconstructor(NULL)
122 ,fClusterIter(NULL)
e3cf3d02 123 ,fExB(0.)
124 ,fVD(0.)
125 ,fT0(0.)
126 ,fS2PRF(0.)
127 ,fDiffL(0.)
128 ,fDiffT(0.)
ae4e8b84 129 ,fClusterIdx(0)
7c3eecb8 130 ,fErrorMsg(0)
3e778975 131 ,fN(0)
e3cf3d02 132 ,fDet(-1)
b25a5e9e 133 ,fPt(0.)
e3cf3d02 134 ,fdX(0.)
135 ,fX0(0.)
136 ,fX(0.)
137 ,fY(0.)
138 ,fZ(0.)
139 ,fS2Y(0.)
140 ,fS2Z(0.)
e3cf3d02 141 ,fChi2(0.)
e4f2f73d 142{
143 //
144 // Copy Constructor performing a deep copy
145 //
e3cf3d02 146 if(this != &ref){
147 ref.Copy(*this);
148 }
29b87567 149 SetBit(kOwner, kFALSE);
f29f13a6 150 SetStandAlone(ref.IsStandAlone());
fbb2ea06 151}
d9950a5a 152
0906e73e 153
e4f2f73d 154//____________________________________________________________________
155AliTRDseedV1& AliTRDseedV1::operator=(const AliTRDseedV1 &ref)
156{
157 //
158 // Assignment Operator using the copy function
159 //
160
29b87567 161 if(this != &ref){
162 ref.Copy(*this);
163 }
221ab7e0 164 SetBit(kOwner, kFALSE);
165
29b87567 166 return *this;
e4f2f73d 167}
168
169//____________________________________________________________________
170AliTRDseedV1::~AliTRDseedV1()
171{
172 //
173 // Destructor. The RecoParam object belongs to the underlying tracker.
174 //
175
29b87567 176 //printf("I-AliTRDseedV1::~AliTRDseedV1() : Owner[%s]\n", IsOwner()?"YES":"NO");
e4f2f73d 177
e3cf3d02 178 if(IsOwner()) {
8d2bec9e 179 for(int itb=0; itb<kNclusters; itb++){
29b87567 180 if(!fClusters[itb]) continue;
181 //AliInfo(Form("deleting c %p @ %d", fClusters[itb], itb));
182 delete fClusters[itb];
4d6aee34 183 fClusters[itb] = NULL;
29b87567 184 }
e3cf3d02 185 }
e4f2f73d 186}
187
188//____________________________________________________________________
189void AliTRDseedV1::Copy(TObject &ref) const
190{
191 //
192 // Copy function
193 //
194
29b87567 195 //AliInfo("");
196 AliTRDseedV1 &target = (AliTRDseedV1 &)ref;
197
4d6aee34 198 target.fkReconstructor = fkReconstructor;
199 target.fClusterIter = NULL;
e3cf3d02 200 target.fExB = fExB;
201 target.fVD = fVD;
202 target.fT0 = fT0;
203 target.fS2PRF = fS2PRF;
204 target.fDiffL = fDiffL;
205 target.fDiffT = fDiffT;
ae4e8b84 206 target.fClusterIdx = 0;
7c3eecb8 207 target.fErrorMsg = fErrorMsg;
3e778975 208 target.fN = fN;
ae4e8b84 209 target.fDet = fDet;
b25a5e9e 210 target.fPt = fPt;
29b87567 211 target.fdX = fdX;
e3cf3d02 212 target.fX0 = fX0;
213 target.fX = fX;
214 target.fY = fY;
215 target.fZ = fZ;
216 target.fS2Y = fS2Y;
217 target.fS2Z = fS2Z;
e3cf3d02 218 target.fChi2 = fChi2;
29b87567 219
8d2bec9e 220 memcpy(target.fIndexes, fIndexes, kNclusters*sizeof(Int_t));
221 memcpy(target.fClusters, fClusters, kNclusters*sizeof(AliTRDcluster*));
2eb10c34 222 memcpy(target.fPad, fPad, 4*sizeof(Float_t));
e3cf3d02 223 target.fYref[0] = fYref[0]; target.fYref[1] = fYref[1];
224 target.fZref[0] = fZref[0]; target.fZref[1] = fZref[1];
225 target.fYfit[0] = fYfit[0]; target.fYfit[1] = fYfit[1];
226 target.fZfit[0] = fZfit[0]; target.fZfit[1] = fZfit[1];
1f97f376 227 memcpy(target.fdEdx, fdEdx, kNdEdxSlices*sizeof(Float_t));
e3cf3d02 228 memcpy(target.fProb, fProb, AliPID::kSPECIES*sizeof(Float_t));
229 memcpy(target.fLabels, fLabels, 3*sizeof(Int_t));
16cca13f 230 memcpy(target.fRefCov, fRefCov, 7*sizeof(Double_t));
68f9b6bd 231 target.fC[0] = fC[0]; target.fC[1] = fC[1];
e3cf3d02 232 memcpy(target.fCov, fCov, 3*sizeof(Double_t));
29b87567 233
e3cf3d02 234 TObject::Copy(ref);
e4f2f73d 235}
236
0906e73e 237
238//____________________________________________________________
4ecadb52 239void AliTRDseedV1::Init(const AliRieman *rieman)
240{
241// Initialize this tracklet using the riemann fit information
242
243
244 fZref[0] = rieman->GetZat(fX0);
245 fZref[1] = rieman->GetDZat(fX0);
246 fYref[0] = rieman->GetYat(fX0);
247 fYref[1] = rieman->GetDYat(fX0);
248 if(fkReconstructor && fkReconstructor->IsHLT()){
249 fRefCov[0] = 1;
250 fRefCov[2] = 10;
251 }else{
252 fRefCov[0] = rieman->GetErrY(fX0);
253 fRefCov[2] = rieman->GetErrZ(fX0);
254 }
255 fC[0] = rieman->GetC();
256 fChi2 = rieman->GetChi2();
257}
258
259
260//____________________________________________________________
33ab3872 261Bool_t AliTRDseedV1::Init(const AliTRDtrackV1 *track)
0906e73e 262{
263// Initialize this tracklet using the track information
264//
265// Parameters:
266// track - the TRD track used to initialize the tracklet
267//
268// Detailed description
269// The function sets the starting point and direction of the
270// tracklet according to the information from the TRD track.
271//
272// Caution
273// The TRD track has to be propagated to the beginning of the
274// chamber where the tracklet will be constructed
275//
276
29b87567 277 Double_t y, z;
278 if(!track->GetProlongation(fX0, y, z)) return kFALSE;
16cca13f 279 Update(track);
29b87567 280 return kTRUE;
0906e73e 281}
282
bcb6fb78 283
e3cf3d02 284//_____________________________________________________________________________
980d5a2a 285void AliTRDseedV1::Reset(Option_t *opt)
e3cf3d02 286{
980d5a2a 287//
288// Reset seed. If option opt="c" is given only cluster arrays are cleared.
289//
290 for(Int_t ic=kNclusters; ic--;) fIndexes[ic] = -1;
291 memset(fClusters, 0, kNclusters*sizeof(AliTRDcluster*));
560e5c05 292 fN=0; SetBit(kRowCross, kFALSE);
980d5a2a 293 if(strcmp(opt, "c")==0) return;
294
e3cf3d02 295 fExB=0.;fVD=0.;fT0=0.;fS2PRF=0.;
296 fDiffL=0.;fDiffT=0.;
3e778975 297 fClusterIdx=0;
7c3eecb8 298 fErrorMsg = 0;
dd8059a8 299 fDet=-1;
b25a5e9e 300 fPt=0.;
e3cf3d02 301 fdX=0.;fX0=0.; fX=0.; fY=0.; fZ=0.;
302 fS2Y=0.; fS2Z=0.;
68f9b6bd 303 fC[0]=0.; fC[1]=0.;
304 fChi2 = 0.;
e3cf3d02 305
2eb10c34 306 memset(fPad, 0, 4*sizeof(Float_t));
e3cf3d02 307 fYref[0] = 0.; fYref[1] = 0.;
308 fZref[0] = 0.; fZref[1] = 0.;
309 fYfit[0] = 0.; fYfit[1] = 0.;
310 fZfit[0] = 0.; fZfit[1] = 0.;
1f97f376 311 memset(fdEdx, 0, kNdEdxSlices*sizeof(Float_t));
e3cf3d02 312 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = -1.;
313 fLabels[0]=-1; fLabels[1]=-1; // most freq MC labels
314 fLabels[2]=0; // number of different labels for tracklet
16cca13f 315 memset(fRefCov, 0, 7*sizeof(Double_t));
e3cf3d02 316 // covariance matrix [diagonal]
317 // default sy = 200um and sz = 2.3 cm
318 fCov[0] = 4.e-4; fCov[1] = 0.; fCov[2] = 5.3;
319}
320
bcb6fb78 321//____________________________________________________________________
16cca13f 322void AliTRDseedV1::Update(const AliTRDtrackV1 *trk)
b1957d3c 323{
324 // update tracklet reference position from the TRD track
b1957d3c 325
e3cf3d02 326 Double_t fSnp = trk->GetSnp();
327 Double_t fTgl = trk->GetTgl();
b25a5e9e 328 fPt = trk->Pt();
bfd20868 329 Double_t norm =1./TMath::Sqrt((1.-fSnp)*(1.+fSnp));
1fd9389f 330 fYref[1] = fSnp*norm;
331 fZref[1] = fTgl*norm;
b1957d3c 332 SetCovRef(trk->GetCovariance());
333
334 Double_t dx = trk->GetX() - fX0;
335 fYref[0] = trk->GetY() - dx*fYref[1];
336 fZref[0] = trk->GetZ() - dx*fZref[1];
337}
338
e3cf3d02 339//_____________________________________________________________________________
340void AliTRDseedV1::UpdateUsed()
341{
342 //
f29f13a6 343 // Calculate number of used clusers in the tracklet
e3cf3d02 344 //
345
3e778975 346 Int_t nused = 0, nshared = 0;
8d2bec9e 347 for (Int_t i = kNclusters; i--; ) {
e3cf3d02 348 if (!fClusters[i]) continue;
3e778975 349 if(fClusters[i]->IsUsed()){
350 nused++;
351 } else if(fClusters[i]->IsShared()){
352 if(IsStandAlone()) nused++;
353 else nshared++;
354 }
e3cf3d02 355 }
3e778975 356 SetNUsed(nused);
357 SetNShared(nshared);
e3cf3d02 358}
359
360//_____________________________________________________________________________
361void AliTRDseedV1::UseClusters()
362{
363 //
364 // Use clusters
365 //
f29f13a6 366 // In stand alone mode:
367 // Clusters which are marked as used or shared from another track are
368 // removed from the tracklet
369 //
370 // In barrel mode:
371 // - Clusters which are used by another track become shared
372 // - Clusters which are attached to a kink track become shared
373 //
e3cf3d02 374 AliTRDcluster **c = &fClusters[0];
8d2bec9e 375 for (Int_t ic=kNclusters; ic--; c++) {
e3cf3d02 376 if(!(*c)) continue;
f29f13a6 377 if(IsStandAlone()){
378 if((*c)->IsShared() || (*c)->IsUsed()){
b82b4de1 379 if((*c)->IsShared()) SetNShared(GetNShared()-1);
380 else SetNUsed(GetNUsed()-1);
4d6aee34 381 (*c) = NULL;
f29f13a6 382 fIndexes[ic] = -1;
3e778975 383 SetN(GetN()-1);
3e778975 384 continue;
f29f13a6 385 }
3e778975 386 } else {
f29f13a6 387 if((*c)->IsUsed() || IsKink()){
3e778975 388 (*c)->SetShared();
389 continue;
f29f13a6 390 }
391 }
392 (*c)->Use();
e3cf3d02 393 }
394}
395
396
f29f13a6 397
b1957d3c 398//____________________________________________________________________
bcb6fb78 399void AliTRDseedV1::CookdEdx(Int_t nslices)
400{
401// Calculates average dE/dx for all slices and store them in the internal array fdEdx.
402//
403// Parameters:
404// nslices : number of slices for which dE/dx should be calculated
405// Output:
406// store results in the internal array fdEdx. This can be accessed with the method
407// AliTRDseedV1::GetdEdx()
408//
409// Detailed description
410// Calculates average dE/dx for all slices. Depending on the PID methode
411// the number of slices can be 3 (LQ) or 8(NN).
3ee48d6e 412// The calculation of dQ/dl are done using the tracklet fit results (see AliTRDseedV1::GetdQdl(Int_t))
bcb6fb78 413//
414// The following effects are included in the calculation:
415// 1. calibration values for t0 and vdrift (using x coordinate to calculate slice)
416// 2. cluster sharing (optional see AliTRDrecoParam::SetClusterSharing())
417// 3. cluster size
418//
419
1f97f376 420 memset(fdEdx, 0, kNdEdxSlices*sizeof(Float_t));
e73abf77 421 const Double_t kDriftLength = (.5 * AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
29b87567 422
0d80a563 423 AliTRDcluster *c(NULL);
424 for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++){
425 if(!(c = fClusters[ic]) && !(c = fClusters[ic+kNtb])) continue;
426 Float_t dx = TMath::Abs(fX0 - c->GetX());
427
29b87567 428 // Filter clusters for dE/dx calculation
0d80a563 429
29b87567 430 // 1.consider calibration effects for slice determination
0d80a563 431 Int_t slice;
432 if(dx<kDriftLength){ // TODO should be replaced by c->IsInChamber()
433 slice = Int_t(dx * nslices / kDriftLength);
434 } else slice = c->GetX() < fX0 ? nslices-1 : 0;
435
e73abf77 436
29b87567 437 // 2. take sharing into account
3e778975 438 Float_t w = /*c->IsShared() ? .5 :*/ 1.;
0d80a563 439
29b87567 440 // 3. take into account large clusters TODO
441 //w *= c->GetNPads() > 3 ? .8 : 1.;
85b917f6 442
0d80a563 443 //CHECK !!!
444 fdEdx[slice] += w * GetdQdl(ic); //fdQdl[ic];
445 } // End of loop over clusters
bcb6fb78 446}
447
e3cf3d02 448//_____________________________________________________________________________
449void AliTRDseedV1::CookLabels()
450{
451 //
452 // Cook 2 labels for seed
453 //
454
455 Int_t labels[200];
456 Int_t out[200];
457 Int_t nlab = 0;
8d2bec9e 458 for (Int_t i = 0; i < kNclusters; i++) {
e3cf3d02 459 if (!fClusters[i]) continue;
460 for (Int_t ilab = 0; ilab < 3; ilab++) {
461 if (fClusters[i]->GetLabel(ilab) >= 0) {
462 labels[nlab] = fClusters[i]->GetLabel(ilab);
463 nlab++;
464 }
465 }
466 }
467
fac58f00 468 fLabels[2] = AliMathBase::Freq(nlab,labels,out,kTRUE);
e3cf3d02 469 fLabels[0] = out[0];
470 if ((fLabels[2] > 1) && (out[3] > 1)) fLabels[1] = out[2];
471}
472
2eb10c34 473//____________________________________________________________
474Float_t AliTRDseedV1::GetAnodeWireOffset(Float_t zt)
475{
4ecadb52 476// Find position inside the amplification cell for reading drift velocity map
477
2eb10c34 478 Float_t d = fPad[3] - zt;
479 if(d<0.){
480 AliError(Form("Fail AnodeWireOffset calculation z0[%+7.2f] zt[%+7.2f] d[%+7.2f].", fPad[3], zt, d));
481 return 0.125;
482 }
483 d -= ((Int_t)(2 * d)) / 2.0;
484 if(d > 0.25) d = 0.5 - d;
485 return d;
486}
487
e3cf3d02 488
d937ad7a 489//____________________________________________________________________
9dcc64cc 490Float_t AliTRDseedV1::GetCharge(Bool_t useOutliers) const
85b917f6 491{
492// Computes total charge attached to tracklet. If "useOutliers" is set clusters
493// which are not in chamber are also used (default false)
494
495 AliTRDcluster *c(NULL); Float_t qt(0.);
496 for(int ic=0; ic<kNclusters; ic++){
497 if(!(c=fClusters[ic])) continue;
06cbae53 498 if(!c->IsInChamber() && !useOutliers) continue;
85b917f6 499 qt += TMath::Abs(c->GetQ());
500 }
501 return qt;
502}
503
504//____________________________________________________________________
15a4c6d0 505Int_t AliTRDseedV1::GetChargeGaps(Float_t sz[kNtb], Float_t pos[kNtb], Int_t isz[kNtb]) const
506{
507// Find number, size and position of charge gaps (consecutive missing time bins).
508// Returns the number of gaps and fills their size in input array "sz" and position in array "pos"
509
510 Bool_t gap(kFALSE);
511 Int_t n(0);
512 Int_t ipos[kNtb]; memset(isz, 0, kNtb*sizeof(Int_t));memset(ipos, 0, kNtb*sizeof(Int_t));
513 for(int ic(0); ic<kNtb; ic++){
514 if(fClusters[ic] || fClusters[ic+kNtb]){
515 if(gap) n++;
516 continue;
517 }
518 gap = kTRUE;
519 isz[n]++;
520 ipos[n] = ic;
521 }
522 if(!n) return 0;
523
524 // write calibrated values
525 AliTRDcluster fake;
526 for(Int_t igap(0); igap<n; igap++){
527 sz[igap] = isz[igap]*fVD/AliTRDCommonParam::Instance()->GetSamplingFrequency();
528 fake.SetPadTime(ipos[igap]);
529 pos[igap] = fake.GetXloc(fT0, fVD);
530 if(isz[igap]>1){
531 fake.SetPadTime(ipos[igap]-isz[igap]+1);
532 pos[igap] += fake.GetXloc(fT0, fVD);
533 pos[igap] /= 2.;
534 }
535 }
536 return n;
537}
538
539
540//____________________________________________________________________
e1bcf0af 541Bool_t AliTRDseedV1::GetEstimatedCrossPoint(Float_t &x, Float_t &z) const
542{
543// Algorithm to estimate cross point in the x-z plane for pad row cross tracklets.
544// Returns true in case of success.
545 if(!IsRowCross()) return kFALSE;
546
547 x=0.; z=0.;
548 AliTRDcluster *c(NULL);
549 // Find radial range for first row
550 Float_t x1[] = {0., 1.e3};
551 for(int ic=0; ic<kNtb; ic++){
552 if(!(c=fClusters[ic])) continue;
553 if(!c->IsInChamber()) continue;
554 if(c->GetX() <= x1[1]) x1[1] = c->GetX();
555 if(c->GetX() >= x1[0]) x1[0] = c->GetX();
556 z=c->GetZ();
557 }
558 if((x1[0] - x1[1])<1.e-5) return kFALSE;
559
560 // Find radial range for second row
561 Bool_t kZ(kFALSE);
562 Float_t x2[] = {0., 1.e3};
563 for(int ic=kNtb; ic<kNclusters; ic++){
564 if(!(c=fClusters[ic])) continue;
565 if(!c->IsInChamber()) continue;
566 if(c->GetX() <= x2[1]) x2[1] = c->GetX();
567 if(c->GetX() >= x2[0]) x2[0] = c->GetX();
568 if(!kZ){
569 z+=c->GetZ();
570 z*=0.5;
571 kZ=kTRUE;
572 }
573 }
574 if((x2[0] - x2[1])<1.e-5) return kFALSE;
575
576 // Find intersection of the 2 radial regions
577 x = 0.5*((x1[0]+x1[1] > x2[0]+x2[1]) ? (x1[1]+x2[0]) : (x1[0]+x2[1]));
578 return kTRUE;
579}
580
581//____________________________________________________________________
a0bb5615 582Float_t AliTRDseedV1::GetQperTB(Int_t tb) const
583{
584 //
585 // Charge of the clusters at timebin
586 //
33ab3872 587 Float_t q = 0;
a0bb5615 588 if(fClusters[tb] /*&& fClusters[tb]->IsInChamber()*/)
33ab3872 589 q += TMath::Abs(fClusters[tb]->GetQ());
a0bb5615 590 if(fClusters[tb+kNtb] /*&& fClusters[tb+kNtb]->IsInChamber()*/)
33ab3872 591 q += TMath::Abs(fClusters[tb+kNtb]->GetQ());
592 return q/TMath::Sqrt(1. + fYref[1]*fYref[1] + fZref[1]*fZref[1]);
a0bb5615 593}
594
595//____________________________________________________________________
9dcc64cc 596Float_t AliTRDseedV1::GetdQdl() const
597{
598// Calculate total charge / tracklet length for 1D PID
599//
600 Float_t Q = GetCharge(kTRUE);
601 return Q/TMath::Sqrt(1. + fYref[1]*fYref[1] + fZref[1]*fZref[1]);
602}
603
604//____________________________________________________________________
0b433f72 605Float_t AliTRDseedV1::GetdQdl(Int_t ic, Float_t *dl) const
bcb6fb78 606{
3ee48d6e 607// Using the linear approximation of the track inside one TRD chamber (TRD tracklet)
608// the charge per unit length can be written as:
609// BEGIN_LATEX
500851ab 610// #frac{dq}{dl} = #frac{q_{c}}{dx * #sqrt{1 + #(){#frac{dy}{dx}}^{2}_{fit} + #(){#frac{dz}{dx}}^{2}_{ref}}}
3ee48d6e 611// END_LATEX
612// where qc is the total charge collected in the current time bin and dx is the length
0b433f72 613// of the time bin.
614// The following correction are applied :
615// - charge : pad row cross corrections
616// [diffusion and TRF assymetry] TODO
617// - dx : anisochronity, track inclination - see Fit and AliTRDcluster::GetXloc()
618// and AliTRDcluster::GetYloc() for the effects taken into account
3ee48d6e 619//
0fa1a8ee 620//Begin_Html
621//<img src="TRD/trackletDQDT.gif">
622//End_Html
623// In the picture the energy loss measured on the tracklet as a function of drift time [left] and respectively
624// drift length [right] for different particle species is displayed.
3ee48d6e 625// Author : Alex Bercuci <A.Bercuci@gsi.de>
626//
627 Float_t dq = 0.;
5d401b45 628 // check whether both clusters are inside the chamber
629 Bool_t hasClusterInChamber = kFALSE;
630 if(fClusters[ic] && fClusters[ic]->IsInChamber()){
631 hasClusterInChamber = kTRUE;
1742f24c 632 dq += TMath::Abs(fClusters[ic]->GetQ());
b30d8c09 633 }
634 if(fClusters[ic+kNtb] && fClusters[ic+kNtb]->IsInChamber()){
5d401b45 635 hasClusterInChamber = kTRUE;
636 dq += TMath::Abs(fClusters[ic+kNtb]->GetQ());
1742f24c 637 }
5d401b45 638 if(!hasClusterInChamber) return 0.;
0b433f72 639 if(dq<1.e-3) return 0.;
3ee48d6e 640
a2abcbc5 641 Double_t dx = fdX;
642 if(ic-1>=0 && ic+1<kNtb){
643 Float_t x2(0.), x1(0.);
5d401b45 644 // try to estimate upper radial position (find the cluster which is inside the chamber)
645 if(fClusters[ic-1] && fClusters[ic-1]->IsInChamber()) x2 = fClusters[ic-1]->GetX();
646 else if(fClusters[ic-1+kNtb] && fClusters[ic-1+kNtb]->IsInChamber()) x2 = fClusters[ic-1+kNtb]->GetX();
647 else if(fClusters[ic] && fClusters[ic]->IsInChamber()) x2 = fClusters[ic]->GetX()+fdX;
a2abcbc5 648 else x2 = fClusters[ic+kNtb]->GetX()+fdX;
5d401b45 649 // try to estimate lower radial position (find the cluster which is inside the chamber)
650 if(fClusters[ic+1] && fClusters[ic+1]->IsInChamber()) x1 = fClusters[ic+1]->GetX();
651 else if(fClusters[ic+1+kNtb] && fClusters[ic+1+kNtb]->IsInChamber()) x1 = fClusters[ic+1+kNtb]->GetX();
652 else if(fClusters[ic] && fClusters[ic]->IsInChamber()) x1 = fClusters[ic]->GetX()-fdX;
a2abcbc5 653 else x1 = fClusters[ic+kNtb]->GetX()-fdX;
654
655 dx = .5*(x2 - x1);
656 }
0b433f72 657 dx *= TMath::Sqrt(1. + fYfit[1]*fYfit[1] + fZref[1]*fZref[1]);
0b433f72 658 if(dl) (*dl) = dx;
283604d2 659 if(dx>1.e-9) return dq/dx;
660 else return 0.;
bcb6fb78 661}
662
0b433f72 663//____________________________________________________________
664Float_t AliTRDseedV1::GetMomentum(Float_t *err) const
665{
666// Returns momentum of the track after update with the current tracklet as:
667// BEGIN_LATEX
668// p=#frac{1}{1/p_{t}} #sqrt{1+tgl^{2}}
669// END_LATEX
670// and optionally the momentum error (if err is not null).
671// The estimated variance of the momentum is given by:
672// BEGIN_LATEX
673// #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})
674// END_LATEX
675// which can be simplified to
676// BEGIN_LATEX
677// #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}
678// END_LATEX
679//
680
681 Double_t p = fPt*TMath::Sqrt(1.+fZref[1]*fZref[1]);
0b433f72 682 if(err){
1f97f376 683 Double_t p2 = p*p;
684 Double_t tgl2 = fZref[1]*fZref[1];
685 Double_t pt2 = fPt*fPt;
686 Double_t s2 =
0b433f72 687 p2*tgl2*pt2*pt2*fRefCov[4]
688 -2.*p2*fZref[1]*fPt*pt2*fRefCov[5]
689 +p2*pt2*fRefCov[6];
690 (*err) = TMath::Sqrt(s2);
691 }
692 return p;
693}
694
9dcc64cc 695
b453ef55 696//____________________________________________________________________
5c5d503a 697Int_t AliTRDseedV1::GetTBoccupancy() const
b453ef55 698{
5c5d503a 699// Returns no. of TB occupied by clusters
b453ef55 700
701 Int_t n(0);
5c5d503a 702 for(int ic(0); ic<kNtb; ic++){
703 if(!fClusters[ic] && !fClusters[ic+kNtb]) continue;
b453ef55 704 n++;
705 }
5c5d503a 706 return n;
707}
b453ef55 708
5c5d503a 709//____________________________________________________________________
710Int_t AliTRDseedV1::GetTBcross() const
711{
712// Returns no. of TB occupied by 2 clusters for pad row cross tracklets
713
714 if(!IsRowCross()) return 0;
715 Int_t n(0);
716 for(int ic(0); ic<kNtb; ic++){
717 if(fClusters[ic] && fClusters[ic+kNtb]) n++;
718 }
719 return n;
b453ef55 720}
0b433f72 721
0906e73e 722//____________________________________________________________________
3e778975 723Float_t* AliTRDseedV1::GetProbability(Bool_t force)
0906e73e 724{
3e778975 725 if(!force) return &fProb[0];
4d6aee34 726 if(!CookPID()) return NULL;
3e778975 727 return &fProb[0];
728}
729
730//____________________________________________________________
731Bool_t AliTRDseedV1::CookPID()
732{
0906e73e 733// Fill probability array for tracklet from the DB.
734//
735// Parameters
736//
737// Output
4d6aee34 738// returns pointer to the probability array and NULL if missing DB access
0906e73e 739//
2a3191bb 740// Retrieve PID probabilities for e+-, mu+-, K+-, pi+- and p+- from the DB according to tracklet information:
741// - estimated momentum at tracklet reference point
742// - dE/dx measurements
743// - tracklet length
744// - TRD layer
745// According to the steering settings specified in the reconstruction one of the following methods are used
746// - Neural Network [default] - option "nn"
747// - 2D Likelihood - option "!nn"
0906e73e 748
1f97f376 749 AliWarning(Form("Obsolete function. Use AliTRDPIDResponse::GetResponse() instead."));
750
0906e73e 751 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
752 if (!calibration) {
753 AliError("No access to calibration data");
3e778975 754 return kFALSE;
0906e73e 755 }
756
4d6aee34 757 if (!fkReconstructor) {
3a039a31 758 AliError("Reconstructor not set.");
3e778975 759 return kFALSE;
4ba1d6ae 760 }
761
0906e73e 762 // Retrieve the CDB container class with the parametric detector response
4d6aee34 763 const AliTRDCalPID *pd = calibration->GetPIDObject(fkReconstructor->GetPIDMethod());
0906e73e 764 if (!pd) {
765 AliError("No access to AliTRDCalPID object");
3e778975 766 return kFALSE;
0906e73e 767 }
10f75631 768
29b87567 769 // calculate tracklet length TO DO
560e5c05 770 Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())/ TMath::Sqrt((1.0 - GetSnp()*GetSnp()) / (1.0 + GetTgl()*GetTgl()));
0906e73e 771
772 //calculate dE/dx
9ded305e 773 CookdEdx(AliTRDCalPID::kNSlicesNN);
774 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 775
0906e73e 776 // Sets the a priori probabilities
11d80e40 777 Bool_t kPIDNN(fkReconstructor->GetPIDMethod()==AliTRDpidUtil::kNN);
f83cd814 778 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++)
11d80e40 779 fProb[ispec] = pd->GetProbability(ispec, GetMomentum(), &fdEdx[0], length, kPIDNN?GetPlane():fkReconstructor->GetRecoParam()->GetPIDLQslices());
f301a656 780
3e778975 781 return kTRUE;
0906e73e 782}
783
784//____________________________________________________________________
e4f2f73d 785Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
786{
787 //
788 // Returns a quality measurement of the current seed
789 //
790
dd8059a8 791 Float_t zcorr = kZcorr ? GetTilt() * (fZfit[0] - fZref[0]) : 0.;
29b87567 792 return
3e778975 793 .5 * TMath::Abs(18.0 - GetN())
29b87567 794 + 10.* TMath::Abs(fYfit[1] - fYref[1])
795 + 5. * TMath::Abs(fYfit[0] - fYref[0] + zcorr)
dd8059a8 796 + 2. * TMath::Abs(fZfit[0] - fZref[0]) / GetPadLength();
e4f2f73d 797}
798
799//____________________________________________________________________
d937ad7a 800void AliTRDseedV1::GetCovAt(Double_t x, Double_t *cov) const
0906e73e 801{
d937ad7a 802// Computes covariance in the y-z plane at radial point x (in tracking coordinates)
803// and returns the results in the preallocated array cov[3] as :
804// cov[0] = Var(y)
805// cov[1] = Cov(yz)
806// cov[2] = Var(z)
807//
808// Details
809//
810// For the linear transformation
811// BEGIN_LATEX
812// Y = T_{x} X^{T}
813// END_LATEX
814// The error propagation has the general form
815// BEGIN_LATEX
816// C_{Y} = T_{x} C_{X} T_{x}^{T}
817// END_LATEX
818// We apply this formula 2 times. First to calculate the covariance of the tracklet
819// at point x we consider:
820// BEGIN_LATEX
821// T_{x} = (1 x); X=(y0 dy/dx); C_{X}=#(){#splitline{Var(y0) Cov(y0, dy/dx)}{Cov(y0, dy/dx) Var(dy/dx)}}
822// END_LATEX
823// and secondly to take into account the tilt angle
824// BEGIN_LATEX
825// T_{#alpha} = #(){#splitline{cos(#alpha) __ sin(#alpha)}{-sin(#alpha) __ cos(#alpha)}}; X=(y z); C_{X}=#(){#splitline{Var(y) 0}{0 Var(z)}}
826// END_LATEX
827//
828// using simple trigonometrics one can write for this last case
829// BEGIN_LATEX
830// 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})}}
831// END_LATEX
832// which can be aproximated for small alphas (2 deg) with
833// BEGIN_LATEX
834// 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}}}
835// END_LATEX
836//
837// before applying the tilt rotation we also apply systematic uncertainties to the tracklet
838// position which can be tunned from outside via the AliTRDrecoParam::SetSysCovMatrix(). They might
839// account for extra misalignment/miscalibration uncertainties.
840//
841// Author :
842// Alex Bercuci <A.Bercuci@gsi.de>
843// Date : Jan 8th 2009
844//
b1957d3c 845
846
d937ad7a 847 Double_t xr = fX0-x;
848 Double_t sy2 = fCov[0] +2.*xr*fCov[1] + xr*xr*fCov[2];
b72f4eaf 849 Double_t sz2 = fS2Z;
850 //GetPadLength()*GetPadLength()/12.;
0906e73e 851
d937ad7a 852 // insert systematic uncertainties
4e6eb3f6 853 /* //RS
4d6aee34 854 if(fkReconstructor){
bb2db46c 855 Double_t sys[15]; memset(sys, 0, 15*sizeof(Double_t));
4d6aee34 856 fkReconstructor->GetRecoParam()->GetSysCovMatrix(sys);
bb2db46c 857 sy2 += sys[0];
858 sz2 += sys[1];
859 }
4e6eb3f6 860 */
2eb10c34 861 // rotate covariance matrix if no RC
862 if(!IsRowCross()){
863 Double_t t2 = GetTilt()*GetTilt();
864 Double_t correction = 1./(1. + t2);
865 cov[0] = (sy2+t2*sz2)*correction;
866 cov[1] = GetTilt()*(sz2 - sy2)*correction;
867 cov[2] = (t2*sy2+sz2)*correction;
868 } else {
a59510a9 869 cov[0] = sy2; cov[1] = 0.; cov[2] = sz2;
2eb10c34 870 }
871
872 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 873}
eb38ed55 874
bb2db46c 875//____________________________________________________________
66765e8e 876Int_t AliTRDseedV1::GetCovSqrt(const Double_t * const c, Double_t *d)
bb2db46c 877{
878// Helper function to calculate the square root of the covariance matrix.
879// The input matrix is stored in the vector c and the result in the vector d.
41b7c7b6 880// Both arrays have to be initialized by the user with at least 3 elements. Return negative in case of failure.
bb2db46c 881//
ec3f0161 882// For calculating the square root of the symmetric matrix c
883// the following relation is used:
bb2db46c 884// BEGIN_LATEX
ec3f0161 885// C^{1/2} = VD^{1/2}V^{-1}
bb2db46c 886// END_LATEX
41b7c7b6 887// with V being the matrix with the n eigenvectors as columns.
ec3f0161 888// In case C is symmetric the followings are true:
889// - matrix D is diagonal with the diagonal given by the eigenvalues of C
41b7c7b6 890// - V = V^{-1}
bb2db46c 891//
892// Author A.Bercuci <A.Bercuci@gsi.de>
893// Date Mar 19 2009
894
66765e8e 895 const Double_t kZero(1.e-20);
4d6aee34 896 Double_t l[2], // eigenvalues
897 v[3]; // eigenvectors
bb2db46c 898 // the secular equation and its solution :
899 // (c[0]-L)(c[2]-L)-c[1]^2 = 0
900 // L^2 - L*Tr(c)+DET(c) = 0
901 // L12 = [Tr(c) +- sqrt(Tr(c)^2-4*DET(c))]/2
4d6aee34 902 Double_t tr = c[0]+c[2], // trace
903 det = c[0]*c[2]-c[1]*c[1]; // determinant
66765e8e 904 if(TMath::Abs(det)<kZero) return 1;
4d6aee34 905 Double_t dd = TMath::Sqrt(tr*tr - 4*det);
66765e8e 906 l[0] = .5*(tr + dd*(c[0]>c[2]?-1.:1.));
907 l[1] = .5*(tr + dd*(c[0]>c[2]?1.:-1.));
908 if(l[0]<kZero || l[1]<kZero) return 2;
41b7c7b6 909 // the sym V matrix
910 // | v00 v10|
911 // | v10 v11|
66765e8e 912 Double_t den = (l[0]-c[0])*(l[0]-c[0])+c[1]*c[1];
913 if(den<kZero){ // almost diagonal
914 v[0] = TMath::Sign(0., c[1]);
915 v[1] = TMath::Sign(1., (l[0]-c[0]));
916 v[2] = TMath::Sign(0., c[1]*(l[0]-c[0])*(l[1]-c[2]));
917 } else {
918 Double_t tmp = 1./TMath::Sqrt(den);
919 v[0] = c[1]* tmp;
920 v[1] = (l[0]-c[0])*tmp;
921 if(TMath::Abs(l[1]-c[2])<kZero) v[2] = TMath::Sign(v[0]*(l[0]-c[0])/kZero, (l[1]-c[2]));
922 else v[2] = v[0]*(l[0]-c[0])/(l[1]-c[2]);
923 }
41b7c7b6 924 // the VD^{1/2}V is:
4d6aee34 925 l[0] = TMath::Sqrt(l[0]); l[1] = TMath::Sqrt(l[1]);
926 d[0] = v[0]*v[0]*l[0]+v[1]*v[1]*l[1];
927 d[1] = v[0]*v[1]*l[0]+v[1]*v[2]*l[1];
928 d[2] = v[1]*v[1]*l[0]+v[2]*v[2]*l[1];
bb2db46c 929
66765e8e 930 return 0;
bb2db46c 931}
932
933//____________________________________________________________
4d6aee34 934Double_t AliTRDseedV1::GetCovInv(const Double_t * const c, Double_t *d)
bb2db46c 935{
936// Helper function to calculate the inverse of the covariance matrix.
937// The input matrix is stored in the vector c and the result in the vector d.
938// Both arrays have to be initialized by the user with at least 3 elements
939// The return value is the determinant or 0 in case of singularity.
940//
941// Author A.Bercuci <A.Bercuci@gsi.de>
942// Date Mar 19 2009
943
4d6aee34 944 Double_t det = c[0]*c[2] - c[1]*c[1];
945 if(TMath::Abs(det)<1.e-20) return 0.;
946 Double_t invDet = 1./det;
947 d[0] = c[2]*invDet;
948 d[1] =-c[1]*invDet;
949 d[2] = c[0]*invDet;
950 return det;
bb2db46c 951}
0906e73e 952
d937ad7a 953//____________________________________________________________________
b72f4eaf 954UShort_t AliTRDseedV1::GetVolumeId() const
955{
4ecadb52 956// Returns geometry volume id by delegation
957
fbe11be7 958 for(Int_t ic(0);ic<kNclusters; ic++){
959 if(fClusters[ic]) return fClusters[ic]->GetVolumeId();
960 }
961 return 0;
b72f4eaf 962}
963
964
965//____________________________________________________________________
e3cf3d02 966void AliTRDseedV1::Calibrate()
d937ad7a 967{
e3cf3d02 968// Retrieve calibration and position parameters from OCDB.
969// The following information are used
d937ad7a 970// - detector index
e3cf3d02 971// - column and row position of first attached cluster. If no clusters are attached
972// to the tracklet a random central chamber position (c=70, r=7) will be used.
973//
974// The following information is cached in the tracklet
975// t0 (trigger delay)
976// drift velocity
977// PRF width
978// omega*tau = tg(a_L)
979// diffusion coefficients (longitudinal and transversal)
d937ad7a 980//
981// Author :
982// Alex Bercuci <A.Bercuci@gsi.de>
983// Date : Jan 8th 2009
984//
eb38ed55 985
d937ad7a 986 AliCDBManager *cdb = AliCDBManager::Instance();
987 if(cdb->GetRun() < 0){
988 AliError("OCDB manager not properly initialized");
989 return;
990 }
0906e73e 991
e3cf3d02 992 AliTRDcalibDB *calib = AliTRDcalibDB::Instance();
993 AliTRDCalROC *vdROC = calib->GetVdriftROC(fDet),
994 *t0ROC = calib->GetT0ROC(fDet);;
995 const AliTRDCalDet *vdDet = calib->GetVdriftDet();
996 const AliTRDCalDet *t0Det = calib->GetT0Det();
d937ad7a 997
998 Int_t col = 70, row = 7;
999 AliTRDcluster **c = &fClusters[0];
3e778975 1000 if(GetN()){
d937ad7a 1001 Int_t ic = 0;
8d2bec9e 1002 while (ic<kNclusters && !(*c)){ic++; c++;}
d937ad7a 1003 if(*c){
1004 col = (*c)->GetPadCol();
1005 row = (*c)->GetPadRow();
1006 }
1007 }
3a039a31 1008
e17f4785 1009 fT0 = (t0Det->GetValue(fDet) + t0ROC->GetValue(col,row)) / AliTRDCommonParam::Instance()->GetSamplingFrequency();
e3cf3d02 1010 fVD = vdDet->GetValue(fDet) * vdROC->GetValue(col, row);
1011 fS2PRF = calib->GetPRFWidth(fDet, col, row); fS2PRF *= fS2PRF;
1012 fExB = AliTRDCommonParam::Instance()->GetOmegaTau(fVD);
1013 AliTRDCommonParam::Instance()->GetDiffCoeff(fDiffL,
1014 fDiffT, fVD);
903326c1 1015 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));
1016
1017
e3cf3d02 1018 SetBit(kCalib, kTRUE);
0906e73e 1019}
1020
0906e73e 1021//____________________________________________________________________
29b87567 1022void AliTRDseedV1::SetOwner()
0906e73e 1023{
29b87567 1024 //AliInfo(Form("own [%s] fOwner[%s]", own?"YES":"NO", fOwner?"YES":"NO"));
1025
1026 if(TestBit(kOwner)) return;
8d2bec9e 1027 for(int ic=0; ic<kNclusters; ic++){
29b87567 1028 if(!fClusters[ic]) continue;
1029 fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
1030 }
1031 SetBit(kOwner);
0906e73e 1032}
1033
eb2b4f91 1034//____________________________________________________________
4ecadb52 1035void AliTRDseedV1::SetPadPlane(AliTRDpadPlane * const p)
eb2b4f91 1036{
1037// Shortcut method to initialize pad geometry.
2eb10c34 1038 fPad[0] = p->GetLengthIPad();
1039 fPad[1] = p->GetWidthIPad();
1040 fPad[2] = TMath::Tan(TMath::DegToRad()*p->GetTiltingAngle());
1041 fPad[3] = p->GetRow0() + p->GetAnodeWireOffset();
eb2b4f91 1042}
1043
1044
9dcc64cc 1045
e4f2f73d 1046//____________________________________________________________________
9dcc64cc 1047Bool_t AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *const chamber, Bool_t tilt, Bool_t chgPos, Int_t ev)
e4f2f73d 1048{
1fd9389f 1049//
1050// Projective algorithm to attach clusters to seeding tracklets. The following steps are performed :
1051// 1. Collapse x coordinate for the full detector plane
1052// 2. truncated mean on y (r-phi) direction
1053// 3. purge clusters
1054// 4. truncated mean on z direction
1055// 5. purge clusters
1056//
1057// Parameters
1058// - chamber : pointer to tracking chamber container used to search the tracklet
1059// - tilt : switch for tilt correction during road building [default true]
803dc399 1060// - chgPos : mark same[kFALSE] and opposite[kTRUE] sign tracks with respect to Bz field sign [default true]
1061// - ev : event number for debug purposes [default = -1]
1fd9389f 1062// Output
1063// - true : if tracklet found successfully. Failure can happend because of the following:
1064// -
1065// Detailed description
9dcc64cc 1066//
1fd9389f 1067// We start up by defining the track direction in the xy plane and roads. The roads are calculated based
8a7ff53c 1068// on tracking information (variance in the r-phi direction) and estimated variance of the standard
1069// clusters (see AliTRDcluster::SetSigmaY2()) corrected for tilt (see GetCovAt()). From this the road is
1070// BEGIN_LATEX
500851ab 1071// 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 1072// r_{z} = 1.5*L_{pad}
1073// END_LATEX
1fd9389f 1074//
4b755889 1075// Author : Alexandru Bercuci <A.Bercuci@gsi.de>
9dcc64cc 1076// Debug : level = 2 for calibration
1077// level = 3 for visualization in the track SR
1078// level = 4 for full visualization including digit level
1fd9389f 1079
fc0882f3 1080 const AliTRDrecoParam* const recoParam = fkReconstructor->GetRecoParam(); //the dynamic cast in GetRecoParam is slow, so caching the pointer to it
1081
1082 if(!recoParam){
560e5c05 1083 AliError("Tracklets can not be used without a valid RecoParam.");
29b87567 1084 return kFALSE;
1085 }
9dcc64cc 1086 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1087 if (!calibration) {
1088 AliError("No access to calibration data");
1089 return kFALSE;
1090 }
1091 // Retrieve the CDB container class with the parametric likelihood
1092 const AliTRDCalTrkAttach *attach = calibration->GetAttachObject();
1093 if (!attach) {
1094 AliError("No usable AttachClusters calib object.");
1095 return kFALSE;
1096 }
1097
b1957d3c 1098 // Initialize reco params for this tracklet
1099 // 1. first time bin in the drift region
a2abcbc5 1100 Int_t t0 = 14;
fc0882f3 1101 Int_t kClmin = Int_t(recoParam->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
9dcc64cc 1102 Int_t kTBmin = 4;
29b87567 1103
9dcc64cc 1104 Double_t sysCov[5]; recoParam->GetSysCovMatrix(sysCov);
8a7ff53c 1105 Double_t s2yTrk= fRefCov[0],
1106 s2yCl = 0.,
1107 s2zCl = GetPadLength()*GetPadLength()/12.,
1108 syRef = TMath::Sqrt(s2yTrk),
1109 t2 = GetTilt()*GetTilt();
29b87567 1110 //define roads
9dcc64cc 1111 const Double_t kroady = 3.; //recoParam->GetRoad1y();
1112 const Double_t kroadz = GetPadLength() * recoParam->GetRoadzMultiplicator() + 1.;
8a7ff53c 1113 // define probing cluster (the perfect cluster) and default calibration
1114 Short_t sig[] = {0, 0, 10, 30, 10, 0,0};
1115 AliTRDcluster cp(fDet, 6, 75, 0, sig, 0);
560e5c05 1116 if(fkReconstructor->IsHLT()) cp.SetRPhiMethod(AliTRDcluster::kCOG);
1117 if(!IsCalibrated()) Calibrate();
8a7ff53c 1118
803dc399 1119/* Int_t kroadyShift(0);
9dcc64cc 1120 Float_t bz(AliTrackerBase::GetBz());
1121 if(TMath::Abs(bz)>2.){
1122 if(bz<0.) kroadyShift = chgPos ? +1 : -1;
1123 else kroadyShift = chgPos ? -1 : +1;
803dc399 1124 }*/
1125 AliDebug(4, Form("\n syTrk[cm]=%4.2f dydxTrk[deg]=%+6.2f Chg[%c] rY[cm]=%4.2f rZ[cm]=%5.2f TC[%c]", syRef, TMath::ATan(fYref[1])*TMath::RadToDeg(), chgPos?'+':'-', kroady, kroadz, tilt?'y':'n'));
9dcc64cc 1126 Double_t phiTrk(TMath::ATan(fYref[1])),
1127 thtTrk(TMath::ATan(fZref[1]));
29b87567 1128
1129 // working variables
b1957d3c 1130 const Int_t kNrows = 16;
4b755889 1131 const Int_t kNcls = 3*kNclusters; // buffer size
9dcc64cc 1132 TObjArray clst[kNrows];
3044dfe5 1133 Bool_t blst[kNrows][kNcls];
9dcc64cc 1134 Double_t cond[4],
1135 dx, dy, dz,
1136 yt, zt,
1137 zc[kNrows],
1138 xres[kNrows][kNcls], yres[kNrows][kNcls], zres[kNrows][kNcls], s2y[kNrows][kNcls];
4b755889 1139 Int_t idxs[kNrows][kNcls], ncl[kNrows], ncls = 0;
b1957d3c 1140 memset(ncl, 0, kNrows*sizeof(Int_t));
9dcc64cc 1141 memset(zc, 0, kNrows*sizeof(Double_t));
1142 memset(idxs, 0, kNrows*kNcls*sizeof(Int_t));
1143 memset(xres, 0, kNrows*kNcls*sizeof(Double_t));
4b755889 1144 memset(yres, 0, kNrows*kNcls*sizeof(Double_t));
9dcc64cc 1145 memset(zres, 0, kNrows*kNcls*sizeof(Double_t));
1146 memset(s2y, 0, kNrows*kNcls*sizeof(Double_t));
3044dfe5 1147 memset(blst, 0, kNrows*kNcls*sizeof(Bool_t)); //this is 8 times faster to memset than "memset(clst, 0, kNrows*kNcls*sizeof(AliTRDcluster*))"
b1957d3c 1148
1f97f376 1149 Double_t roady(0.), s2Mean(0.); Int_t ns2Mean(0);
9dcc64cc 1150
1151 // Do cluster projection and pick up cluster candidates
1152 AliTRDcluster *c(NULL);
1153 AliTRDchamberTimeBin *layer(NULL);
b1957d3c 1154 Bool_t kBUFFER = kFALSE;
4b755889 1155 for (Int_t it = 0; it < kNtb; it++) {
b1957d3c 1156 if(!(layer = chamber->GetTB(it))) continue;
29b87567 1157 if(!Int_t(*layer)) continue;
8a7ff53c 1158 // get track projection at layers position
b1957d3c 1159 dx = fX0 - layer->GetX();
1160 yt = fYref[0] - fYref[1] * dx;
1161 zt = fZref[0] - fZref[1] * dx;
9dcc64cc 1162 // get standard cluster error corrected for tilt if selected
8a7ff53c 1163 cp.SetLocalTimeBin(it);
1164 cp.SetSigmaY2(0.02, fDiffT, fExB, dx, -1./*zt*/, fYref[1]);
9dcc64cc 1165 s2yCl = cp.GetSigmaY2() + sysCov[0]; if(!tilt) s2yCl = (s2yCl + t2*s2zCl)/(1.+t2);
1166 if(TMath::Abs(it-12)<7){ s2Mean += cp.GetSigmaY2(); ns2Mean++;}
1167 // get estimated road in r-phi direction
1168 roady = TMath::Min(3.*TMath::Sqrt(12.*(s2yTrk + s2yCl)), kroady);
1169
1170 AliDebug(5, Form("\n"
1171 " %2d xd[cm]=%6.3f yt[cm]=%7.2f zt[cm]=%8.2f\n"
1172 " syTrk[um]=%6.2f syCl[um]=%6.2f syClTlt[um]=%6.2f\n"
1173 " Ry[mm]=%f"
1174 , it, dx, yt, zt
1175 , 1.e4*TMath::Sqrt(s2yTrk), 1.e4*TMath::Sqrt(cp.GetSigmaY2()+sysCov[0]), 1.e4*TMath::Sqrt(s2yCl)
1176 , 1.e1*roady));
1177
1178 // get clusters from layer
1179 cond[0] = yt/*+0.5*kroadyShift*kroady*/; cond[2] = roady;
b1957d3c 1180 cond[1] = zt; cond[3] = kroadz;
9dcc64cc 1181 Int_t n=0, idx[6]; layer->GetClusters(cond, idx, n, 6);
b1957d3c 1182 for(Int_t ic = n; ic--;){
1183 c = (*layer)[idx[ic]];
9dcc64cc 1184 dx = fX0 - c->GetX();
1185 yt = fYref[0] - fYref[1] * dx;
1186 zt = fZref[0] - fZref[1] * dx;
1187 dz = zt - c->GetZ();
1188 dy = yt - (c->GetY() + (tilt ? (GetTilt() * dz) : 0.));
b1957d3c 1189 Int_t r = c->GetPadRow();
9dcc64cc 1190 clst[r].AddAtAndExpand(c, ncl[r]);
3044dfe5 1191 blst[r][ncl[r]] = kTRUE;
b1957d3c 1192 idxs[r][ncl[r]] = idx[ic];
9dcc64cc 1193 zres[r][ncl[r]] = dz/GetPadLength();
b1957d3c 1194 yres[r][ncl[r]] = dy;
9dcc64cc 1195 xres[r][ncl[r]] = dx;
1196 zc[r] = c->GetZ();
1197 // TODO temporary solution to avoid divercences in error parametrization
1198 s2y[r][ncl[r]] = TMath::Min(c->GetSigmaY2()+sysCov[0], 0.025);
1199 AliDebug(5, Form(" -> dy[cm]=%+7.4f yc[cm]=%7.2f row[%d] idx[%2d]", dy, c->GetY(), r, ncl[r]));
b1957d3c 1200 ncl[r]++; ncls++;
1201
4b755889 1202 if(ncl[r] >= kNcls) {
560e5c05 1203 AliWarning(Form("Cluster candidates row[%d] reached buffer limit[%d]. Some may be lost.", r, kNcls));
b1957d3c 1204 kBUFFER = kTRUE;
29b87567 1205 break;
1206 }
1207 }
b1957d3c 1208 if(kBUFFER) break;
29b87567 1209 }
ee8fb199 1210 if(ncls<kClmin){
560e5c05 1211 AliDebug(1, Form("CLUSTERS FOUND %d LESS THAN THRESHOLD %d.", ncls, kClmin));
7c3eecb8 1212 SetErrorMsg(kAttachClFound);
9dcc64cc 1213 for(Int_t ir(kNrows);ir--;) clst[ir].Clear();
ee8fb199 1214 return kFALSE;
1215 }
9dcc64cc 1216 if(ns2Mean<kTBmin){
1217 AliDebug(1, Form("CLUSTERS IN TimeBins %d LESS THAN THRESHOLD %d.", ns2Mean, kTBmin));
1218 SetErrorMsg(kAttachClFound);
1219 for(Int_t ir(kNrows);ir--;) clst[ir].Clear();
1220 return kFALSE;
1221 }
1f97f376 1222 s2Mean /= ns2Mean; //sMean = TMath::Sqrt(s2Mean);
9dcc64cc 1223 //Double_t sRef(TMath::Sqrt(s2Mean+s2yTrk)); // reference error parameterization
1224
1225 // organize row candidates
1226 Int_t idxRow[kNrows], nrc(0); Double_t zresRow[kNrows];
1227 for(Int_t ir(0); ir<kNrows; ir++){
1228 idxRow[ir]=-1; zresRow[ir] = 999.;
1229 if(!ncl[ir]) continue;
1230 // get mean z resolution
1231 dz = 0.; for(Int_t ic = ncl[ir]; ic--;) dz += zres[ir][ic]; dz/=ncl[ir];
1232 // insert row
1233 idxRow[nrc] = ir; zresRow[nrc] = TMath::Abs(dz); nrc++;
1234 }
1235 AliDebug(4, Form("Found %d clusters in %d rows. Sorting ...", ncls, nrc));
1236
1237 // sort row candidates
1238 if(nrc>=2){
1239 if(nrc==2){
1240 if(zresRow[0]>zresRow[1]){ // swap
1241 Int_t itmp=idxRow[1]; idxRow[1] = idxRow[0]; idxRow[0] = itmp;
1242 Double_t dtmp=zresRow[1]; zresRow[1] = zresRow[0]; zresRow[0] = dtmp;
1243 }
1244 if(TMath::Abs(idxRow[1] - idxRow[0]) != 1){
1245 SetErrorMsg(kAttachRowGap);
1246 AliDebug(2, Form("Rows attached not continuous. Select first candidate.\n"
1247 " row[%2d] Ncl[%2d] <dz>[cm]=%+8.2f row[%2d] Ncl[%2d] <dz>[cm]=%+8.2f",
0a62661e 1248 idxRow[0], ncl[idxRow[0]], zresRow[0], idxRow[1], idxRow[1]<0?0:ncl[idxRow[1]], zresRow[1]));
9dcc64cc 1249 nrc=1; idxRow[1] = -1; zresRow[1] = 999.;
1250 }
1251 } else {
1252 Int_t idx0[kNrows];
1253 TMath::Sort(nrc, zresRow, idx0, kFALSE);
1254 nrc = 3; // select only maximum first 3 candidates
1255 Int_t iatmp[] = {-1, -1, -1}; Double_t datmp[] = {999., 999., 999.};
1256 for(Int_t irc(0); irc<nrc; irc++){
1257 iatmp[irc] = idxRow[idx0[irc]];
1258 datmp[irc] = zresRow[idx0[irc]];
1259 }
1260 idxRow[0] = iatmp[0]; zresRow[0] = datmp[0];
1261 idxRow[1] = iatmp[1]; zresRow[1] = datmp[1];
1262 idxRow[2] = iatmp[2]; zresRow[2] = datmp[2]; // temporary
1263 if(TMath::Abs(idxRow[1] - idxRow[0]) != 1){
1264 SetErrorMsg(kAttachRowGap);
1265 AliDebug(2, Form("Rows attached not continuous. Turn on selection.\n"
1266 "row[%2d] Ncl[%2d] <dz>[cm]=%+8.2f\n"
1267 "row[%2d] Ncl[%2d] <dz>[cm]=%+8.2f\n"
1268 "row[%2d] Ncl[%2d] <dz>[cm]=%+8.2f",
1269 idxRow[0], ncl[idxRow[0]], zresRow[0],
1270 idxRow[1], ncl[idxRow[1]], zresRow[1],
1271 idxRow[2], ncl[idxRow[2]], zresRow[2]));
1272 if(TMath::Abs(idxRow[0] - idxRow[2]) == 1){ // select second candidate
1273 AliDebug(2, "Solved ! Remove second candidate.");
1274 nrc = 2;
1275 idxRow[1] = idxRow[2]; zresRow[1] = zresRow[2]; // swap
1276 idxRow[2] = -1; zresRow[2] = 999.; // remove
1277 } else if(TMath::Abs(idxRow[1] - idxRow[2]) == 1){
1278 if(ncl[idxRow[1]]+ncl[idxRow[2]] > ncl[idxRow[0]]){
1279 AliDebug(2, "Solved ! Remove first candidate.");
1280 nrc = 2;
1281 idxRow[0] = idxRow[1]; zresRow[0] = zresRow[1]; // swap
1282 idxRow[1] = idxRow[2]; zresRow[1] = zresRow[2]; // swap
1283 } else {
1284 AliDebug(2, "Solved ! Remove second and third candidate.");
1285 nrc = 1;
1286 idxRow[1] = -1; zresRow[1] = 999.; // remove
1287 idxRow[2] = -1; zresRow[2] = 999.; // remove
1288 }
1289 } else {
1290 AliDebug(2, "Unsolved !!! Remove second and third candidate.");
1291 nrc = 1;
1292 idxRow[1] = -1; zresRow[1] = 999.; // remove
1293 idxRow[2] = -1; zresRow[2] = 999.; // remove
1294 }
1295 } else { // remove temporary candidate
1296 nrc = 2;
1297 idxRow[2] = -1; zresRow[2] = 999.;
b1957d3c 1298 }
b1957d3c 1299 }
29b87567 1300 }
9dcc64cc 1301 AliDebug(4, Form("Sorted row candidates:\n"
1302 " row[%2d] Ncl[%2d] <dz>[cm]=%+8.2f row[%2d] Ncl[%2d] <dz>[cm]=%+8.2f"
dbb2e0a7 1303 , idxRow[0], ncl[idxRow[0]], zresRow[0], idxRow[1], idxRow[1]<0?0:ncl[idxRow[1]], zresRow[1]));
9dcc64cc 1304
1305 // initialize debug streamer
1306 TTreeSRedirector *pstreamer(NULL);
121f4b48 1307 if(recoParam->GetStreamLevel(AliTRDrecoParam::kTracker) > 3 && fkReconstructor->IsDebugStreaming()) pstreamer = fkReconstructor->GetDebugStream(AliTRDrecoParam::kTracker);
9dcc64cc 1308 if(pstreamer){
1309 // save config. for calibration
1310 TVectorD vdy[2], vdx[2], vs2[2];
1311 for(Int_t jr(0); jr<nrc; jr++){
1312 Int_t ir(idxRow[jr]);
1313 vdx[jr].ResizeTo(ncl[ir]); vdy[jr].ResizeTo(ncl[ir]); vs2[jr].ResizeTo(ncl[ir]);
1314 for(Int_t ic(ncl[ir]); ic--;){
1315 vdx[jr](ic) = xres[ir][ic];
1316 vdy[jr](ic) = yres[ir][ic];
1317 vs2[jr](ic) = s2y[ir][ic];
1318 }
1319 }
1320 (*pstreamer) << "AttachClusters4"
1321 << "r0=" << idxRow[0]
1322 << "dz0=" << zresRow[0]
1323 << "dx0=" << &vdx[0]
560e5c05 1324 << "dy0=" << &vdy[0]
9dcc64cc 1325 << "s20=" << &vs2[0]
1326 << "r1=" << idxRow[1]
1327 << "dz1=" << zresRow[1]
1328 << "dx1=" << &vdx[1]
560e5c05 1329 << "dy1=" << &vdy[1]
9dcc64cc 1330 << "s21=" << &vs2[1]
560e5c05 1331 << "\n";
9dcc64cc 1332 vdx[0].Clear(); vdy[0].Clear(); vs2[0].Clear();
1333 vdx[1].Clear(); vdy[1].Clear(); vs2[1].Clear();
121f4b48 1334 if(recoParam->GetStreamLevel(AliTRDrecoParam::kTracker) > 4){
2f4384e6 1335 Int_t idx(idxRow[1]);
1336 if(idx<0){
9dcc64cc 1337 for(Int_t ir(0); ir<kNrows; ir++){
1338 if(clst[ir].GetEntries()>0) continue;
1339 idx = ir;
1340 break;
1341 }
2f4384e6 1342 }
9dcc64cc 1343 (*pstreamer) << "AttachClusters5"
1344 << "c0.=" << &clst[idxRow[0]]
1345 << "c1.=" << &clst[idx]
1346 << "\n";
1347 }
560e5c05 1348 }
1349
9dcc64cc 1350//=======================================================================================
1351 // Analyse cluster topology
1352 Double_t f[kNcls], // likelihood factors for segments
1353 r[2][kNcls], // d(dydx) of tracklet candidate with respect to track
1354 xm[2][kNcls], // mean <x>
1355 ym[2][kNcls], // mean <y>
1356 sm[2][kNcls], // mean <s_y>
1357 s[2][kNcls], // sigma_y
2f4384e6 1358 p[2][kNcls], // prob of Gauss
1359 q[2][kNcls]; // charge/segment
9dcc64cc 1360 memset(f, 0, kNcls*sizeof(Double_t));
1361 Int_t index[2][kNcls], n[2][kNcls];
1362 memset(n, 0, 2*kNcls*sizeof(Int_t));
1363 Int_t mts(0), nts[2] = {0, 0}; // no of tracklet segments in row
1364 AliTRDpadPlane *pp(AliTRDtransform::Geometry().GetPadPlane(fDet));
1365 AliTRDtrackletOflHelper helper;
1366 Int_t lyDet(AliTRDgeometry::GetLayer(fDet));
1367 for(Int_t jr(0), n0(0); jr<nrc; jr++){
1368 Int_t ir(idxRow[jr]);
1369 // cluster segmentation
1370 Bool_t kInit(kFALSE);
1371 if(jr==0){
1372 n0 = helper.Init(pp, &clst[ir]); kInit = kTRUE;
1373 if(!n0 || (helper.ClassifyTopology() == AliTRDtrackletOflHelper::kNormal)){
1374 nts[jr] = 1; memset(index[jr], 0, ncl[ir]*sizeof(Int_t));
1375 n[jr][0] = ncl[ir];
560e5c05 1376 }
9dcc64cc 1377 }
1378 if(!n[jr][0]){
1379 nts[jr] = AliTRDtrackletOflHelper::Segmentation(ncl[ir], xres[ir], yres[ir], index[jr]);
1380 for(Int_t ic(ncl[ir]);ic--;) n[jr][index[jr][ic]]++;
1381 }
1382 mts += nts[jr];
1383
1384 // tracklet segment processing
1385 for(Int_t its(0); its<nts[jr]; its++){
1386 if(n[jr][its]<=2) { // don't touch small segments
1387 xm[jr][its] = 0.;ym[jr][its] = 0.;sm[jr][its] = 0.;
1388 for(Int_t ic(ncl[ir]); ic--;){
1389 if(its != index[jr][ic]) continue;
1390 ym[jr][its] += yres[ir][ic];
1391 xm[jr][its] += xres[ir][ic];
1392 sm[jr][its] += TMath::Sqrt(s2y[ir][ic]);
1393 }
1394 if(n[jr][its]==2){ xm[jr][its] *= 0.5; ym[jr][its] *= 0.5; sm[jr][its] *= 0.5;}
1395 xm[jr][its]= fX0 - xm[jr][its];
1396 r[jr][its] = 0.;
1397 s[jr][its] = 1.e-5;
1398 p[jr][its] = 1.;
2f4384e6 1399 q[jr][its] = -1.;
9dcc64cc 1400 continue;
560e5c05 1401 }
9dcc64cc 1402
1403 // for longer tracklet segments
1404 if(!kInit) n0 = helper.Init(pp, &clst[ir], index[jr], its);
2f4384e6 1405 Int_t n1 = helper.GetRMS(r[jr][its], ym[jr][its], s[jr][its], fX0/*xm[jr][its]*/);
1406 p[jr][its] = Double_t(n1)/n0;
9dcc64cc 1407 sm[jr][its] = helper.GetSyMean();
2f4384e6 1408 q[jr][its] = helper.GetQ()/TMath::Sqrt(1. + fYref[1]*fYref[1] + fZref[1]*fZref[1]);
1409 xm[jr][its] = fX0;
9dcc64cc 1410 Double_t dxm= fX0 - xm[jr][its];
2f4384e6 1411 yt = fYref[0] - fYref[1]*dxm;
9dcc64cc 1412 zt = fZref[0] - fZref[1]*dxm;
1413 // correct tracklet fit for tilt
1414 ym[jr][its]+= GetTilt()*(zt - zc[ir]);
1415 r[jr][its] += GetTilt() * fZref[1];
1416 // correct tracklet fit for track position/inclination
2f4384e6 1417 ym[jr][its] = yt - ym[jr][its];
1418 r[jr][its] = (r[jr][its] - fYref[1])/(1+r[jr][its]*fYref[1]);
9dcc64cc 1419 // report inclination in radians
1420 r[jr][its] = TMath::ATan(r[jr][its]);
1421 if(jr) continue; // calculate only for first row likelihoods
1422
1423 f[its] = attach->CookLikelihood(chgPos, lyDet, fPt, phiTrk, n[jr][its], ym[jr][its]/*sRef*/, r[jr][its]*TMath::RadToDeg(), s[jr][its]/sm[jr][its]);
29b87567 1424 }
1425 }
9dcc64cc 1426 AliDebug(4, Form(" Tracklet candidates: row[%2d] = %2d row[%2d] = %2d:", idxRow[0], nts[0], idxRow[1], nts[1]));
1427 if(AliLog::GetDebugLevel("TRD", "AliTRDseedV1")>3){
1428 for(Int_t jr(0); jr<nrc; jr++){
1429 Int_t ir(idxRow[jr]);
1430 for(Int_t its(0); its<nts[jr]; its++){
1431 printf(" segId[%2d] row[%2d] Ncl[%2d] x[cm]=%7.2f dz[pu]=%4.2f dy[mm]=%+7.3f r[deg]=%+6.2f p[%%]=%6.2f s[um]=%7.2f\n",
1432 its, ir, n[jr][its], xm[jr][its], zresRow[jr], 1.e1*ym[jr][its], r[jr][its]*TMath::RadToDeg(), 100.*p[jr][its], 1.e4*s[jr][its]);
1433 }
560e5c05 1434 }
ee8fb199 1435 }
9dcc64cc 1436 if(!pstreamer && recoParam->GetStreamLevel(AliTRDrecoParam::kTracker) > 2 && fkReconstructor->IsDebugStreaming()) pstreamer = fkReconstructor->GetDebugStream(AliTRDrecoParam::kTracker);
1437 if(pstreamer){
1438 // save config. for calibration
1439 TVectorD vidx, vn, vx, vy, vr, vs, vsm, vp, vf;
1440 vidx.ResizeTo(ncl[idxRow[0]]+(idxRow[1]<0?0:ncl[idxRow[1]]));
1441 vn.ResizeTo(mts);
1442 vx.ResizeTo(mts);
1443 vy.ResizeTo(mts);
1444 vr.ResizeTo(mts);
1445 vs.ResizeTo(mts);
1446 vsm.ResizeTo(mts);
1447 vp.ResizeTo(mts);
1448 vf.ResizeTo(mts);
1449 for(Int_t jr(0), jts(0), jc(0); jr<nrc; jr++){
1450 Int_t ir(idxRow[jr]);
1451 for(Int_t its(0); its<nts[jr]; its++, jts++){
1452 vn[jts] = n[jr][its];
1453 vx[jts] = xm[jr][its];
1454 vy[jts] = ym[jr][its];
1455 vr[jts] = r[jr][its];
1456 vs[jts] = s[jr][its];
1457 vsm[jts]= sm[jr][its];
1458 vp[jts] = p[jr][its];
1459 vf[jts] = jr?-1.:f[its];
6ad5e6b2 1460 }
9dcc64cc 1461 for(Int_t ic(0); ic<ncl[ir]; ic++, jc++) vidx[jc] = index[jr][ic];
1462 }
1463 (*pstreamer) << "AttachClusters3"
1464 << "idx=" << &vidx
1465 << "n=" << &vn
1466 << "x=" << &vx
1467 << "y=" << &vy
1468 << "r=" << &vr
1469 << "s=" << &vs
1470 << "sm=" << &vsm
1471 << "p=" << &vp
1472 << "f=" << &vf
1473 << "\n";
1474 }
6ad5e6b2 1475
9dcc64cc 1476//=========================================================
1477 // Get seed tracklet segment
1478 Int_t idx2[kNcls]; memset(idx2, 0, kNcls*sizeof(Int_t)); // seeding indexing
1479 if(nts[0]>1) TMath::Sort(nts[0], f, idx2);
1480 Int_t is(idx2[0]); // seed index
1481 Int_t idxTrklt[kNcls],
1482 kts(0),
1483 nTrklt(n[0][is]);
1484 Double_t fTrklt(f[is]),
1485 rTrklt(r[0][is]),
1486 yTrklt(ym[0][is]),
1487 sTrklt(s[0][is]),
1488 smTrklt(sm[0][is]),
1489 xTrklt(xm[0][is]),
2f4384e6 1490 pTrklt(p[0][is]),
1491 qTrklt(q[0][is]);
9dcc64cc 1492 memset(idxTrklt, 0, kNcls*sizeof(Int_t));
1493 // check seed idx2[0] exit if not found
1494 if(f[is]<1.e-2){
1495 AliDebug(1, Form("Seed seg[%d] row[%2d] n[%2d] f[%f]<0.01.", is, idxRow[0], n[0][is], f[is]));
1496 SetErrorMsg(kAttachClAttach);
1497 if(!pstreamer && recoParam->GetStreamLevel(AliTRDrecoParam::kTracker) > 1 && fkReconstructor->IsDebugStreaming()) pstreamer = fkReconstructor->GetDebugStream(AliTRDrecoParam::kTracker);
1498 if(pstreamer){
1499 UChar_t stat(0);
1500 if(IsKink()) SETBIT(stat, 1);
1501 if(IsStandAlone()) SETBIT(stat, 2);
1502 if(IsRowCross()) SETBIT(stat, 3);
1503 SETBIT(stat, 4); // set error bit
1504 TVectorD vidx; vidx.ResizeTo(1); vidx[0] = is;
1505 (*pstreamer) << "AttachClusters2"
1506 << "stat=" << stat
1507 << "ev=" << ev
1508 << "chg=" << chgPos
1509 << "det=" << fDet
1510 << "x0=" << fX0
1511 << "y0=" << fYref[0]
1512 << "z0=" << fZref[0]
1513 << "phi=" << phiTrk
1514 << "tht=" << thtTrk
1515 << "pt=" << fPt
1516 << "s2Trk=" << s2yTrk
1517 << "s2Cl=" << s2Mean
1518 << "idx=" << &vidx
1519 << "n=" << nTrklt
1520 << "f=" << fTrklt
1521 << "x=" << xTrklt
1522 << "y=" << yTrklt
1523 << "r=" << rTrklt
1524 << "s=" << sTrklt
1525 << "sm=" << smTrklt
1526 << "p=" << pTrklt
2f4384e6 1527 << "q=" << qTrklt
9dcc64cc 1528 << "\n";
1529 }
1530 return kFALSE;
1531 }
2f4384e6 1532 AliDebug(2, Form("Seed seg[%d] row[%2d] n[%2d] dy[%f] r[%+5.2f] s[%+5.2f] f[%5.3f] q[%6.2f]", is, idxRow[0], n[0][is], ym[0][is], r[0][is]*TMath::RadToDeg(), s[0][is]/sm[0][is], f[is], q[0][is]));
9dcc64cc 1533
1534 // save seeding segment in the helper
1535 idxTrklt[kts++] = is;
1536 helper.Init(pp, &clst[idxRow[0]], index[0], is);
1537 AliTRDtrackletOflHelper test; // helper to test segment expantion
1538 Float_t rcLikelihood(0.); SetBit(kRowCross, kFALSE);
1539 Double_t dyRez[kNcls]; Int_t idx3[kNcls];
1540
1541 //=========================================================
1542 // Define filter parameters from OCDB
1543 Int_t kNSgmDy[2]; attach->GetNsgmDy(kNSgmDy[0], kNSgmDy[1]);
1544 Float_t kLikeMinRelDecrease[2]; attach->GetLikeMinRelDecrease(kLikeMinRelDecrease[0], kLikeMinRelDecrease[1]);
1545 Float_t kRClikeLimit(attach->GetRClikeLimit());
1546
1547 //=========================================================
1548 // Try attaching next segments from first row (if any)
1549 if(nts[0]>1){
1550 Int_t jr(0), ir(idxRow[jr]);
1551 // organize secondary sgms. in decreasing order of their distance from seed
1552 memset(dyRez, 0, nts[jr]*sizeof(Double_t));
1553 for(Int_t jts(1); jts<nts[jr]; jts++) {
1554 Int_t its(idx2[jts]);
1555 Double_t rot(TMath::Tan(r[0][is]));
1556 dyRez[its] = TMath::Abs(ym[0][is] - ym[jr][its] + rot*(xm[0][is]-xm[jr][its]));
1557 }
1558 TMath::Sort(nts[jr], dyRez, idx3, kFALSE);
1559 for (Int_t jts(1); jts<nts[jr]; jts++) {
1560 Int_t its(idx3[jts]);
1561 if(dyRez[its] > kNSgmDy[jr]*smTrklt){
1562 AliDebug(2, Form("Reject seg[%d] row[%2d] n[%2d] dy[%f] > %d*s[%f].", its, idxRow[jr], n[jr][its], dyRez[its], kNSgmDy[jr], kNSgmDy[jr]*smTrklt));
1563 continue;
1564 }
1565
1566 test = helper;
1567 Int_t n0 = test.Expand(&clst[ir], index[jr], its);
2f4384e6 1568 Double_t rt, dyt, st, xt, smt, pt, qt, ft;
1569 Int_t n1 = test.GetRMS(rt, dyt, st, fX0/*xt*/);
9dcc64cc 1570 pt = Double_t(n1)/n0;
1571 smt = test.GetSyMean();
2f4384e6 1572 qt = test.GetQ()/TMath::Sqrt(1. + fYref[1]*fYref[1] + fZref[1]*fZref[1]);
1573 xt = fX0;
9dcc64cc 1574 // correct position
1575 Double_t dxm= fX0 - xt;
1576 yt = fYref[0] - fYref[1]*dxm;
1577 zt = fZref[0] - fZref[1]*dxm;
1578 // correct tracklet fit for tilt
1579 dyt+= GetTilt()*(zt - zc[idxRow[0]]);
1580 rt += GetTilt() * fZref[1];
1581 // correct tracklet fit for track position/inclination
2f4384e6 1582 dyt = yt - dyt;
1583 rt = (rt - fYref[1])/(1+rt*fYref[1]);
9dcc64cc 1584 // report inclination in radians
1585 rt = TMath::ATan(rt);
1586
1587 ft = (n0>=2) ? attach->CookLikelihood(chgPos, lyDet, fPt, phiTrk, n0, dyt/*sRef*/, rt*TMath::RadToDeg(), st/smt) : 0.;
1588 Bool_t kAccept(ft>=fTrklt*(1.-kLikeMinRelDecrease[jr]));
1589
1590 AliDebug(2, Form("%s seg[%d] row[%2d] n[%2d] dy[%f] r[%+5.2f] s[%+5.2f] f[%f] < %4.2f*F[%f].",
1591 (kAccept?"Adding":"Reject"), its, idxRow[jr], n0, dyt, rt*TMath::RadToDeg(), st/smt, ft, 1.-kLikeMinRelDecrease[jr], fTrklt*(1.-kLikeMinRelDecrease[jr])));
1592 if(kAccept){
1593 idxTrklt[kts++] = its;
1594 nTrklt = n0;
1595 fTrklt = ft;
1596 rTrklt = rt;
1597 yTrklt = dyt;
1598 sTrklt = st;
1599 smTrklt= smt;
1600 xTrklt = xt;
1601 pTrklt = pt;
2f4384e6 1602 qTrklt = qt;
9dcc64cc 1603 helper.Expand(&clst[ir], index[jr], its);
1604 }
b1957d3c 1605 }
560e5c05 1606 }
9dcc64cc 1607
1608 //=========================================================
1609 // Try attaching next segments from second row (if any)
1610 if(nts[1] && (rcLikelihood = zresRow[0]/zresRow[1]) > kRClikeLimit){
1611 // organize secondaries in decreasing order of their distance from seed
1612 Int_t jr(1), ir(idxRow[jr]);
1613 memset(dyRez, 0, nts[jr]*sizeof(Double_t));
1614 Double_t rot(TMath::Tan(r[0][is]));
1615 for(Int_t jts(0); jts<nts[jr]; jts++) {
1616 dyRez[jts] = TMath::Abs(ym[0][is] - ym[jr][jts] + rot*(xm[0][is]-xm[jr][jts]));
1617 }
1618 TMath::Sort(nts[jr], dyRez, idx3, kFALSE);
1619 for (Int_t jts(0); jts<nts[jr]; jts++) {
1620 Int_t its(idx3[jts]);
1621 if(dyRez[its] > kNSgmDy[jr]*smTrklt){
1622 AliDebug(2, Form("Reject seg[%d] row[%2d] n[%2d] dy[%f] > %d*s[%f].", its, idxRow[jr], n[jr][its], dyRez[its], kNSgmDy[jr], kNSgmDy[jr]*smTrklt));
1623 continue;
1624 }
1625
1626 test = helper;
1627 Int_t n0 = test.Expand(&clst[ir], index[jr], its);
2f4384e6 1628 Double_t rt, dyt, st, xt, smt, pt, qt, ft;
1629 Int_t n1 = test.GetRMS(rt, dyt, st, fX0/*xt*/);
9dcc64cc 1630 pt = Double_t(n1)/n0;
1631 smt = test.GetSyMean();
2f4384e6 1632 qt = test.GetQ()/TMath::Sqrt(1. + fYref[1]*fYref[1] + fZref[1]*fZref[1]);
1633 xt = fX0;
9dcc64cc 1634 // correct position
1635 Double_t dxm= fX0 - xt;
1636 yt = fYref[0] - fYref[1]*dxm;
1637 zt = fZref[0] - fZref[1]*dxm;
1638 // correct tracklet fit for tilt
1639 dyt+= GetTilt()*(zt - zc[idxRow[0]]);
1640 rt += GetTilt() * fZref[1];
1641 // correct tracklet fit for track position/inclination
2f4384e6 1642 dyt = yt - dyt;
1643 rt = (rt - fYref[1])/(1+rt*fYref[1]);
9dcc64cc 1644 // report inclination in radians
1645 rt = TMath::ATan(rt);
1646
1647 ft = (n0>=2) ? attach->CookLikelihood(chgPos, lyDet, fPt, phiTrk, n0, dyt/*sRef*/, rt*TMath::RadToDeg(), st/smt) : 0.;
1648 Bool_t kAccept(ft>=fTrklt*(1.-kLikeMinRelDecrease[jr]));
1649
1650 AliDebug(2, Form("%s seg[%d] row[%2d] n[%2d] dy[%f] r[%+5.2f] s[%+5.2f] f[%f] < %4.2f*F[%f].",
1651 (kAccept?"Adding":"Reject"), its, idxRow[jr], n0, dyt, rt*TMath::RadToDeg(), st/smt, ft, 1.-kLikeMinRelDecrease[jr], fTrklt*(1.-kLikeMinRelDecrease[jr])));
1652 if(kAccept){
1653 idxTrklt[kts++] = its;
1654 nTrklt = n0;
1655 fTrklt = ft;
1656 rTrklt = rt;
1657 yTrklt = dyt;
1658 sTrklt = st;
1659 smTrklt= smt;
1660 xTrklt = xt;
1661 pTrklt = pt;
2f4384e6 1662 qTrklt = qt;
9dcc64cc 1663 helper.Expand(&clst[ir], index[jr], its);
1664 SetBit(kRowCross, kTRUE); // mark pad row crossing
1665 }
1666 }
1667 }
1668 // clear local copy of clusters
1669 for(Int_t ir(0); ir<kNrows; ir++) clst[ir].Clear();
1670
1671 if(!pstreamer && recoParam->GetStreamLevel(AliTRDrecoParam::kTracker) > 1 && fkReconstructor->IsDebugStreaming()) pstreamer = fkReconstructor->GetDebugStream(AliTRDrecoParam::kTracker);
1672 if(pstreamer){
1673 UChar_t stat(0);
1674 if(IsKink()) SETBIT(stat, 1);
1675 if(IsStandAlone()) SETBIT(stat, 2);
1676 if(IsRowCross()) SETBIT(stat, 3);
1677 TVectorD vidx; vidx.ResizeTo(kts);
1678 for(Int_t its(0); its<kts; its++) vidx[its] = idxTrklt[its];
1679 (*pstreamer) << "AttachClusters2"
1680 << "stat=" << stat
1681 << "ev=" << ev
1682 << "chg=" << chgPos
1683 << "det=" << fDet
1684 << "x0=" << fX0
1685 << "y0=" << fYref[0]
1686 << "z0=" << fZref[0]
1687 << "phi=" << phiTrk
1688 << "tht=" << thtTrk
1689 << "pt=" << fPt
1690 << "s2Trk=" << s2yTrk
1691 << "s2Cl=" << s2Mean
1692 << "idx=" << &vidx
1693 << "n=" << nTrklt
2f4384e6 1694 << "q=" << qTrklt
9dcc64cc 1695 << "f=" << fTrklt
1696 << "x=" << xTrklt
1697 << "y=" << yTrklt
1698 << "r=" << rTrklt
1699 << "s=" << sTrklt
1700 << "sm=" << smTrklt
1701 << "p=" << pTrklt
1702 << "\n";
1703 }
1704
1705
1706 //=========================================================
1707 // Store clusters
1708 Int_t nselected(0), nc(0);
1709 TObjArray *selected(helper.GetClusters());
1710 if(!selected || !(nselected = selected->GetEntriesFast())){
1711 AliError("Cluster candidates missing !!!");
1712 SetErrorMsg(kAttachClAttach);
1713 return kFALSE;
1714 }
1715 for(Int_t ic(0); ic<nselected; ic++){
1716 if(!(c = (AliTRDcluster*)selected->At(ic))) continue;
1717 Int_t it(c->GetPadTime()),
1718 jr(Int_t(helper.GetRow() != c->GetPadRow())),
1719 idx(it+kNtb*jr);
1720 if(fClusters[idx]){
1721 AliDebug(1, Form("Multiple clusters/tb for D[%03d] Tb[%02d] Row[%2d]", fDet, it, c->GetPadRow()));
1722 continue; // already booked
1723 }
1724 // TODO proper indexing of clusters !!
1725 fIndexes[idx] = chamber->GetTB(it)->GetGlobalIndex(idxs[idxRow[jr]][ic]);
1726 fClusters[idx] = c;
1727 nc++;
1728 }
1729 AliDebug(2, Form("Clusters Found[%2d] Attached[%2d] RC[%c]", nselected, nc, IsRowCross()?'y':'n'));
b1957d3c 1730
29b87567 1731 // number of minimum numbers of clusters expected for the tracklet
9dcc64cc 1732 if (nc < kClmin){
1733 AliDebug(1, Form("NOT ENOUGH CLUSTERS %d ATTACHED TO THE TRACKLET [min %d] FROM FOUND %d.", nc, kClmin, ncls));
7c3eecb8 1734 SetErrorMsg(kAttachClAttach);
e4f2f73d 1735 return kFALSE;
1736 }
9dcc64cc 1737 SetN(nc);
0906e73e 1738
e3cf3d02 1739 // Load calibration parameters for this tracklet
9dcc64cc 1740 //Calibrate();
b1957d3c 1741
1742 // calculate dx for time bins in the drift region (calibration aware)
a2abcbc5 1743 Float_t x[2] = {0.,0.}; Int_t tb[2]={0,0};
1744 for (Int_t it = t0, irp=0; irp<2 && it < AliTRDtrackerV1::GetNTimeBins(); it++) {
b1957d3c 1745 if(!fClusters[it]) continue;
1746 x[irp] = fClusters[it]->GetX();
a2abcbc5 1747 tb[irp] = fClusters[it]->GetLocalTimeBin();
b1957d3c 1748 irp++;
e3cf3d02 1749 }
d86ed84c 1750 Int_t dtb = tb[1] - tb[0];
1751 fdX = dtb ? (x[0] - x[1]) / dtb : 0.15;
29b87567 1752 return kTRUE;
e4f2f73d 1753}
1754
03cef9b2 1755//____________________________________________________________
1756void AliTRDseedV1::Bootstrap(const AliTRDReconstructor *rec)
1757{
1758// Fill in all derived information. It has to be called after recovery from file or HLT.
1759// The primitive data are
1760// - list of clusters
1761// - detector (as the detector will be removed from clusters)
1762// - position of anode wire (fX0) - temporary
1763// - track reference position and direction
1764// - momentum of the track
1765// - time bin length [cm]
1766//
1767// A.Bercuci <A.Bercuci@gsi.de> Oct 30th 2008
1768//
4d6aee34 1769 fkReconstructor = rec;
03cef9b2 1770 AliTRDgeometry g;
2eb10c34 1771 SetPadPlane(g.GetPadPlane(fDet));
1772
e3cf3d02 1773 //fSnp = fYref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);
1774 //fTgl = fZref[1];
3e778975 1775 Int_t n = 0, nshare = 0, nused = 0;
03cef9b2 1776 AliTRDcluster **cit = &fClusters[0];
8d2bec9e 1777 for(Int_t ic = kNclusters; ic--; cit++){
03cef9b2 1778 if(!(*cit)) return;
3e778975 1779 n++;
1780 if((*cit)->IsShared()) nshare++;
1781 if((*cit)->IsUsed()) nused++;
03cef9b2 1782 }
3e778975 1783 SetN(n); SetNUsed(nused); SetNShared(nshare);
e3cf3d02 1784 Fit();
03cef9b2 1785 CookLabels();
1786 GetProbability();
1787}
1788
1789
e4f2f73d 1790//____________________________________________________________________
2eb10c34 1791Bool_t AliTRDseedV1::Fit(UChar_t opt)
e4f2f73d 1792{
16cca13f 1793//
1794// Linear fit of the clusters attached to the tracklet
1795//
1796// Parameters :
2eb10c34 1797// - opt : switch for tilt pad correction of cluster y position. Options are
1798// 0 no correction [default]
1799// 1 full tilt correction [dz/dx and z0]
1800// 2 pseudo tilt correction [dz/dx from pad-chamber geometry]
1801//
16cca13f 1802// Output :
1803// True if successful
1804//
1805// Detailed description
1806//
1807// Fit in the xy plane
1808//
1fd9389f 1809// The fit is performed to estimate the y position of the tracklet and the track
1810// angle in the bending plane. The clusters are represented in the chamber coordinate
1811// system (with respect to the anode wire - see AliTRDtrackerV1::FollowBackProlongation()
1812// on how this is set). The x and y position of the cluster and also their variances
1813// are known from clusterizer level (see AliTRDcluster::GetXloc(), AliTRDcluster::GetYloc(),
1814// AliTRDcluster::GetSX() and AliTRDcluster::GetSY()).
1815// If gaussian approximation is used to calculate y coordinate of the cluster the position
1816// is recalculated taking into account the track angle. The general formula to calculate the
1817// error of cluster position in the gaussian approximation taking into account diffusion and track
1818// inclination is given for TRD by:
1819// BEGIN_LATEX
1820// #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}
1821// END_LATEX
1822//
1823// Since errors are calculated only in the y directions, radial errors (x direction) are mapped to y
1824// by projection i.e.
1825// BEGIN_LATEX
1826// #sigma_{x|y} = tg(#phi) #sigma_{x}
1827// END_LATEX
1828// and also by the lorentz angle correction
1829//
1830// Fit in the xz plane
1831//
1832// The "fit" is performed to estimate the radial position (x direction) where pad row cross happens.
1833// If no pad row crossing the z position is taken from geometry and radial position is taken from the xy
1834// fit (see below).
1835//
1836// There are two methods to estimate the radial position of the pad row cross:
1837// 1. leading cluster radial position : Here the lower part of the tracklet is considered and the last
1838// cluster registered (at radial x0) on this segment is chosen to mark the pad row crossing. The error
1839// of the z estimate is given by :
1840// BEGIN_LATEX
1841// #sigma_{z} = tg(#theta) #Delta x_{x_{0}}/12
1842// END_LATEX
1843// The systematic errors for this estimation are generated by the following sources:
1844// - no charge sharing between pad rows is considered (sharp cross)
1845// - missing cluster at row cross (noise peak-up, under-threshold signal etc.).
1846//
1847// 2. charge fit over the crossing point : Here the full energy deposit along the tracklet is considered
1848// to estimate the position of the crossing by a fit in the qx plane. The errors in the q directions are
1849// parameterized as s_q = q^2. The systematic errors for this estimation are generated by the following sources:
1850// - no general model for the qx dependence
1851// - physical fluctuations of the charge deposit
1852// - gain calibration dependence
1853//
1854// Estimation of the radial position of the tracklet
16cca13f 1855//
1fd9389f 1856// For pad row cross the radial position is taken from the xz fit (see above). Otherwise it is taken as the
1857// interpolation point of the tracklet i.e. the point where the error in y of the fit is minimum. The error
1858// in the y direction of the tracklet is (see AliTRDseedV1::GetCovAt()):
1859// BEGIN_LATEX
1860// #sigma_{y} = #sigma^{2}_{y_{0}} + 2xcov(y_{0}, dy/dx) + #sigma^{2}_{dy/dx}
1861// END_LATEX
1862// and thus the radial position is:
1863// BEGIN_LATEX
1864// x = - cov(y_{0}, dy/dx)/#sigma^{2}_{dy/dx}
1865// END_LATEX
1866//
1867// Estimation of tracklet position error
1868//
1869// The error in y direction is the error of the linear fit at the radial position of the tracklet while in the z
1870// direction is given by the cluster error or pad row cross error. In case of no pad row cross this is given by:
1871// BEGIN_LATEX
1872// #sigma_{y} = #sigma^{2}_{y_{0}} - 2cov^{2}(y_{0}, dy/dx)/#sigma^{2}_{dy/dx} + #sigma^{2}_{dy/dx}
1873// #sigma_{z} = Pad_{length}/12
1874// END_LATEX
1875// For pad row cross the full error is calculated at the radial position of the crossing (see above) and the error
1876// in z by the width of the crossing region - being a matter of parameterization.
1877// BEGIN_LATEX
1878// #sigma_{z} = tg(#theta) #Delta x_{x_{0}}/12
1879// END_LATEX
1880// In case of no tilt correction (default in the barrel tracking) the tilt is taken into account by the rotation of
1881// the covariance matrix. See AliTRDseedV1::GetCovAt() for details.
1882//
1883// Author
1884// A.Bercuci <A.Bercuci@gsi.de>
e4f2f73d 1885
a723055f 1886 if(!fkReconstructor){
1887 AliError("The tracklet needs the reconstruction setup. Please initialize by SetReconstructor().");
1888 return kFALSE;
1889 }
b72f4eaf 1890 if(!IsCalibrated()) Calibrate();
2eb10c34 1891 if(opt>2){
7e5954f0 1892 AliWarning(Form("Option [%d] outside range [0, 2]. Using default",opt));
2eb10c34 1893 opt=0;
1894 }
e3cf3d02 1895
29b87567 1896 const Int_t kClmin = 8;
2eb10c34 1897 const Float_t kScalePulls = 10.; // factor to scale y pulls - NOT UNDERSTOOD
2f7d6ac8 1898 // get track direction
1899 Double_t y0 = fYref[0];
1900 Double_t dydx = fYref[1];
1901 Double_t z0 = fZref[0];
1902 Double_t dzdx = fZref[1];
ae4e8b84 1903
5f1ae1e7 1904 AliTRDtrackerV1::AliTRDLeastSquare fitterY;
1905 AliTRDtrackerV1::AliTRDLeastSquare fitterZ;
f301a656 1906
29b87567 1907 // book cluster information
8d2bec9e 1908 Double_t qc[kNclusters], xc[kNclusters], yc[kNclusters], zc[kNclusters], sy[kNclusters];
e3cf3d02 1909
2eb10c34 1910 Bool_t tilt(opt==1) // full tilt correction
1911 ,pseudo(opt==2) // pseudo tilt correction
1912 ,rc(IsRowCross()) // row cross candidate
1913 ,kDZDX(IsPrimary());// switch dzdx calculation for barrel primary tracks
1914 Int_t n(0); // clusters used in fit
1915 AliTRDcluster *c(NULL), *cc(NULL), **jc = &fClusters[0];
fc0882f3 1916 const AliTRDrecoParam* const recoParam = fkReconstructor->GetRecoParam(); //the dynamic cast in GetRecoParam is slow, so caching the pointer to it
2eb10c34 1917
1918 const Char_t *tcName[]={"NONE", "FULL", "HALF"};
1919 AliDebug(2, Form("Options : TC[%s] dzdx[%c]", tcName[opt], kDZDX?'Y':'N'));
1920
9dcc64cc 1921
2eb10c34 1922 for (Int_t ic=0; ic<kNclusters; ic++, ++jc) {
1923 xc[ic] = -1.; yc[ic] = 999.; zc[ic] = 999.; sy[ic] = 0.;
9eb2d46c 1924 if(!(c = (*jc))) continue;
29b87567 1925 if(!c->IsInChamber()) continue;
2eb10c34 1926 // compute pseudo tilt correction
1927 if(kDZDX){
1928 fZfit[0] = c->GetZ();
1929 if(rc){
1930 for(Int_t kc=AliTRDseedV1::kNtb; kc<AliTRDseedV1::kNclusters; kc++){
1931 if(!(cc=fClusters[kc])) continue;
1932 if(!cc->IsInChamber()) continue;
1933 fZfit[0] += cc->GetZ(); fZfit[0] *= 0.5;
1934 break;
1935 }
1936 }
1937 fZfit[1] = fZfit[0]/fX0;
1938 if(rc){
1939 fZfit[0] += fZfit[1]*0.5*AliTRDgeometry::CdrHght();
1940 fZfit[1] = fZfit[0]/fX0;
1941 }
1942 kDZDX=kFALSE;
1943 }
9462866a 1944
1f97f376 1945// TODO use this information to adjust cluster error parameterization
1946// Float_t w = 1.;
1947// if(c->GetNPads()>4) w = .5;
1948// if(c->GetNPads()>5) w = .2;
010d62b0 1949
1fd9389f 1950 // cluster charge
dd8059a8 1951 qc[n] = TMath::Abs(c->GetQ());
1fd9389f 1952 // pad row of leading
1953
b72f4eaf 1954 xc[n] = fX0 - c->GetX();
1955
1fd9389f 1956 // Recalculate cluster error based on tracking information
2eb10c34 1957 c->SetSigmaY2(fS2PRF, fDiffT, fExB, xc[n], -1./*zcorr?zt:-1.*/, dydx);
c79857d5 1958 c->SetSigmaZ2(fPad[0]*fPad[0]/12.); // for HLT
1fd9389f 1959 sy[n] = TMath::Sqrt(c->GetSigmaY2());
1960
fc0882f3 1961 yc[n] = recoParam->UseGAUS() ?
1fd9389f 1962 c->GetYloc(y0, sy[n], GetPadWidth()): c->GetY();
1963 zc[n] = c->GetZ();
2eb10c34 1964
1965 //optional r-phi correction
1966 //printf(" n[%2d] yc[%7.5f] ", n, yc[n]);
1967 Float_t correction(0.);
1968 if(tilt) correction = fPad[2]*(xc[n]*dzdx + zc[n] - z0);
1969 else if(pseudo) correction = fPad[2]*(xc[n]*fZfit[1] + zc[n]-fZfit[0]);
1970 yc[n]-=correction;
1971 //printf("corr(%s%s)[%7.5f] yc1[%7.5f]\n", (tilt?"TC":""), (zcorr?"PC":""), correction, yc[n]);
1fd9389f 1972
fbe11be7 1973 AliDebug(5, Form(" tb[%2d] dx[%6.3f] y[%6.2f+-%6.3f]", c->GetLocalTimeBin(), xc[n], yc[n], sy[n]));
903326c1 1974 fitterY.AddPoint(&xc[n], yc[n], sy[n]);
2eb10c34 1975 if(rc) fitterZ.AddPoint(&xc[n], qc[n]*(ic<kNtb?1.:-1.), 1.);
dd8059a8 1976 n++;
29b87567 1977 }
3044dfe5 1978
47d5d320 1979 // to few clusters
c79857d5 1980 if (n < kClmin){
c388cdcb 1981 AliDebug(1, Form("Not enough clusters to fit. Clusters: Attached[%d] Fit[%d].", GetN(), n));
2eb10c34 1982 SetErrorMsg(kFitCl);
c79857d5 1983 return kFALSE;
1984 }
d937ad7a 1985 // fit XY
903326c1 1986 if(!fitterY.Eval()){
c388cdcb 1987 AliDebug(1, "Fit Y failed.");
2eb10c34 1988 SetErrorMsg(kFitFailedY);
903326c1 1989 return kFALSE;
1990 }
5f1ae1e7 1991 fYfit[0] = fitterY.GetFunctionParameter(0);
1992 fYfit[1] = -fitterY.GetFunctionParameter(1);
d937ad7a 1993 // store covariance
5f1ae1e7 1994 Double_t p[3];
1995 fitterY.GetCovarianceMatrix(p);
2eb10c34 1996 fCov[0] = kScalePulls*p[1]; // variance of y0
1997 fCov[1] = kScalePulls*p[2]; // covariance of y0, dydx
1998 fCov[2] = kScalePulls*p[0]; // variance of dydx
b1957d3c 1999 // the ref radial position is set at the minimum of
2000 // the y variance of the tracklet
b72f4eaf 2001 fX = -fCov[1]/fCov[2];
2eb10c34 2002 fS2Y = fCov[0] +2.*fX*fCov[1] + fX*fX*fCov[2];
2003
903326c1 2004 Float_t xs=fX+.5*AliTRDgeometry::CamHght();
2005 if(xs < 0. || xs > AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght()){
2006 AliDebug(1, Form("Ref radial position ouside chamber x[%5.2f].", fX));
2eb10c34 2007 SetErrorMsg(kFitFailedY);
903326c1 2008 return kFALSE;
2009 }
b1957d3c 2010
2eb10c34 2011/* // THE LEADING CLUSTER METHOD for z fit
1fd9389f 2012 Float_t xMin = fX0;
b72f4eaf 2013 Int_t ic=n=kNclusters-1; jc = &fClusters[ic];
1fd9389f 2014 AliTRDcluster *c0 =0x0, **kc = &fClusters[kNtb-1];
2015 for(; ic>kNtb; ic--, --jc, --kc){
2016 if((c0 = (*kc)) && c0->IsInChamber() && (xMin>c0->GetX())) xMin = c0->GetX();
2017 if(!(c = (*jc))) continue;
2018 if(!c->IsInChamber()) continue;
2019 zc[kNclusters-1] = c->GetZ();
2020 fX = fX0 - c->GetX();
2021 }
2022 fZfit[0] = .5*(zc[0]+zc[kNclusters-1]); fZfit[1] = 0.;
2023 // Error parameterization
2024 fS2Z = fdX*fZref[1];
e355f67a 2025 fS2Z *= fS2Z; fS2Z *= 0.2887; // 1/sqrt(12)*/
2026
2eb10c34 2027 // fit QZ
2028 if(opt!=1 && IsRowCross()){
2029 if(!fitterZ.Eval()) SetErrorMsg(kFitFailedZ);
4ecadb52 2030 if(!HasError(kFitFailedZ) && TMath::Abs(fitterZ.GetFunctionParameter(1))>1.e-10){
2eb10c34 2031 // TODO - one has to recalculate xy fit based on
2032 // better knowledge of z position
2033// Double_t x = -fitterZ.GetFunctionParameter(0)/fitterZ.GetFunctionParameter(1);
2034// Double_t z0 = .5*(zc[0]+zc[n-1]);
2035// fZfit[0] = z0 + fZfit[1]*x;
2036// fZfit[1] = fZfit[0]/fX0;
2037// redo fit on xy plane
b72f4eaf 2038 }
c850c351 2039 // temporary external error parameterization
2040 fS2Z = 0.05+0.4*TMath::Abs(fZref[1]); fS2Z *= fS2Z;
2041 // TODO correct formula
2042 //fS2Z = sigma_x*TMath::Abs(fZref[1]);
b1957d3c 2043 } else {
2eb10c34 2044 //fZfit[0] = zc[0] + dzdx*0.5*AliTRDgeometry::CdrHght();
dd8059a8 2045 fS2Z = GetPadLength()*GetPadLength()/12.;
29b87567 2046 }
29b87567 2047 return kTRUE;
e4f2f73d 2048}
2049
e4f2f73d 2050
9dcc64cc 2051//____________________________________________________________________
2052Bool_t AliTRDseedV1::FitRobust(Bool_t chg)
e3cf3d02 2053{
2054//
9dcc64cc 2055// Linear fit of the clusters attached to the tracklet
e3cf3d02 2056//
9dcc64cc 2057// Author
2058// A.Bercuci <A.Bercuci@gsi.de>
e3cf3d02 2059
9dcc64cc 2060 TTreeSRedirector *pstreamer(NULL);
2061 const AliTRDrecoParam* const recoParam = fkReconstructor->GetRecoParam(); if(recoParam->GetStreamLevel(AliTRDrecoParam::kTracker) > 3 && fkReconstructor->IsDebugStreaming()) pstreamer = fkReconstructor->GetDebugStream(AliTRDrecoParam::kTracker);
e3cf3d02 2062
9dcc64cc 2063 // factor to scale y pulls.
2064 // ideally if error parametrization correct this is 1.
2065 //Float_t lyScaler = 1./(AliTRDgeometry::GetLayer(fDet)+1.);
2066 Float_t kScalePulls = 1.;
2067 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
2068 if(!calibration){
2069 AliWarning("No access to calibration data");
2070 } else {
2071 // Retrieve the CDB container class with the parametric likelihood
2072 const AliTRDCalTrkAttach *attach = calibration->GetAttachObject();
2073 if(!attach){
2074 AliWarning("No usable AttachClusters calib object.");
2075 } else {
4e6eb3f6 2076 // kScalePulls = attach->GetScaleCov();//*lyScaler;
9dcc64cc 2077 }
803dc399 2078 // Retrieve chamber status
2079 SetChmbGood(calibration->IsChamberGood(fDet));
2080 if(!IsChmbGood()) kScalePulls*=10.;
e3cf3d02 2081 }
9dcc64cc 2082 Double_t xc[kNclusters], yc[kNclusters], sy[kNclusters];
2083 Int_t n(0), // clusters used in fit
2084 row[]={-1, 0}; // pad row spanned by the tracklet
2085 AliTRDcluster *c(NULL), **jc = &fClusters[0];
2086 for(Int_t ic=0; ic<kNtb; ic++, ++jc) {
2087 if(!(c = (*jc))) continue;
2088 if(!c->IsInChamber()) continue;
2089 if(row[0]<0){
2090 fZfit[0] = c->GetZ();
2091 fZfit[1] = 0.;
2092 row[0] = c->GetPadRow();
e3cf3d02 2093 }
2f4384e6 2094 xc[n] = c->GetX();
9dcc64cc 2095 yc[n] = c->GetY();
2096 sy[n] = c->GetSigmaY2()>0?(TMath::Min(TMath::Sqrt(c->GetSigmaY2()), 0.08)):0.08;
2097 n++;
e3cf3d02 2098 }
9dcc64cc 2099 Double_t corr = fPad[2]*fPad[0];
e3cf3d02 2100
9dcc64cc 2101 for(Int_t ic=kNtb; ic<kNclusters; ic++, ++jc) {
2102 if(!(c = (*jc))) continue;
2103 if(!c->IsInChamber()) continue;
2104 if(row[1]==0) row[1] = c->GetPadRow() - row[0];
2f4384e6 2105 xc[n] = c->GetX();
9dcc64cc 2106 yc[n] = c->GetY() + corr*row[1];
2107 sy[n] = c->GetSigmaY2()>0?(TMath::Min(TMath::Sqrt(c->GetSigmaY2()), 0.08)):0.08;
2108 n++;
e3cf3d02 2109 }
9dcc64cc 2110 UChar_t status(0);
2f4384e6 2111 Double_t par[3] = {0.,0.,fX0}, cov[3];
9dcc64cc 2112 if(!AliTRDtrackletOflHelper::Fit(n, xc, yc, sy, par, 1.5, cov)){
2113 AliDebug(1, Form("Tracklet fit failed D[%03d].", fDet));
2114 SetErrorMsg(kFitCl);
2115 return kFALSE;
e3cf3d02 2116 }
9dcc64cc 2117 fYfit[0] = par[0];
2f4384e6 2118 fYfit[1] = par[1];
9dcc64cc 2119 // store covariance
2120 fCov[0] = kScalePulls*cov[0]; // variance of y0
2121 fCov[1] = kScalePulls*cov[2]; // covariance of y0, dydx
2122 fCov[2] = kScalePulls*cov[1]; // variance of dydx
2123 // the ref radial position is set at the minimum of
2124 // the y variance of the tracklet
2f4384e6 2125 fX = 0.;//-fCov[1]/fCov[2];
9dcc64cc 2126 // check radial position
2127 Float_t xs=fX+.5*AliTRDgeometry::CamHght();
2128 if(xs < 0. || xs > AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght()){
2129 AliDebug(1, Form("Ref radial position x[%5.2f] ouside D[%3d].", fX, fDet));
2130 SetErrorMsg(kFitFailedY);
2131 return kFALSE;
e3cf3d02 2132 }
9dcc64cc 2133 fS2Y = fCov[0] + fX*fCov[1];
2134 fS2Z = fPad[0]*fPad[0]/12.;
2135 AliDebug(2, Form("[I] x[cm]=%6.2f y[cm]=%+5.2f z[cm]=%+6.2f dydx[deg]=%+5.2f sy[um]=%6.2f sz[cm]=%6.2f", GetX(), GetY(), GetZ(), TMath::ATan(fYfit[1])*TMath::RadToDeg(), TMath::Sqrt(fS2Y)*1.e4, TMath::Sqrt(fS2Z)));
2136 if(IsRowCross()){
2137 Float_t x,z;
2138 if(!GetEstimatedCrossPoint(x,z)){
2f4384e6 2139 AliDebug(2, Form("Failed(I) getting crossing point D[%03d].", fDet));
9dcc64cc 2140 SetErrorMsg(kFitFailedY);
2141 return kTRUE;
2142 }
2f4384e6 2143 //if(IsPrimary()){
2144 fZfit[0] = fX0*z/x;
9dcc64cc 2145 fZfit[1] = z/x;
2146 fS2Z = 0.05+0.4*TMath::Abs(fZfit[1]); fS2Z *= fS2Z;
2f4384e6 2147 //}
2148 AliDebug(2, Form("s2y[%f] s2z[%f]", fS2Y, fS2Z));
9dcc64cc 2149 AliDebug(2, Form("[II] x[cm]=%6.2f y[cm]=%+5.2f z[cm]=%+6.2f dydx[deg]=%+5.2f sy[um]=%6.2f sz[um]=%6.2f dzdx[deg]=%+5.2f", GetX(), GetY(), GetZ(), TMath::ATan(fYfit[1])*TMath::RadToDeg(), TMath::Sqrt(fS2Y)*1.e4, TMath::Sqrt(fS2Z)*1.e4, TMath::ATan(fZfit[1])*TMath::RadToDeg()));
e3cf3d02 2150 }
e3cf3d02 2151
9dcc64cc 2152 if(pstreamer){
2153 Float_t x= fX0 -fX,
2154 y = GetY(),
2155 yt = fYref[0]-fX*fYref[1];
2156 SETBIT(status, 2);
2157 TVectorD vcov(3); vcov[0]=cov[0];vcov[1]=cov[1];vcov[2]=cov[2];
2158 Double_t sm(0.), chi2(0.), tmp, dy[kNclusters];
2159 for(Int_t ic(0); ic<n; ic++){
2160 sm += sy[ic];
2f4384e6 2161 dy[ic] = yc[ic]-(fYfit[0]+(xc[ic]-fX0)*fYfit[1]); tmp = dy[ic]/sy[ic];
9dcc64cc 2162 chi2 += tmp*tmp;
2163 }
2164 sm /= n; chi2 = TMath::Sqrt(chi2);
2165 Double_t m(0.), s(0.);
2166 AliMathBase::EvaluateUni(n, dy, m, s, 0);
2167 (*pstreamer) << "FitRobust4"
2168 << "stat=" << status
2169 << "chg=" << chg
2170 << "ncl=" << n
2171 << "det=" << fDet
2172 << "x0=" << fX0
2173 << "y0=" << fYfit[0]
2174 << "x=" << x
2175 << "y=" << y
2176 << "dydx=" << fYfit[1]
2177 << "pt=" << fPt
2178 << "yt=" << yt
2179 << "dydxt="<< fYref[1]
2180 << "cov=" << &vcov
2181 << "chi2=" << chi2
2182 << "sm=" << sm
2183 << "ss=" << s
2184 << "\n";
2185 }
2186 return kTRUE;
2187}
e3cf3d02 2188
e4f2f73d 2189//___________________________________________________________________
203967fc 2190void AliTRDseedV1::Print(Option_t *o) const
e4f2f73d 2191{
2192 //
2193 // Printing the seedstatus
2194 //
2195
b72f4eaf 2196 AliInfo(Form("Det[%3d] X0[%7.2f] Pad{L[%5.2f] W[%5.2f] Tilt[%+6.2f]}", fDet, fX0, GetPadLength(), GetPadWidth(), GetTilt()));
dd8059a8 2197 AliInfo(Form("N[%2d] Nused[%2d] Nshared[%2d] [%d]", GetN(), GetNUsed(), GetNShared(), fN));
b72f4eaf 2198 AliInfo(Form("FLAGS : RC[%c] Kink[%c] SA[%c]", IsRowCross()?'y':'n', IsKink()?'y':'n', IsStandAlone()?'y':'n'));
525f399b 2199 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 2200
2201 Double_t cov[3], x=GetX();
2202 GetCovAt(x, cov);
2203 AliInfo(" | x[cm] | y[cm] | z[cm] | dydx | dzdx |");
2204 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]));
21aad919 2205 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 2206 AliInfo(Form("P / Pt [GeV/c] = %f / %f", GetMomentum(), fPt));
68f9b6bd 2207 if(IsStandAlone()) AliInfo(Form("C Rieman / Vertex [1/cm] = %f / %f", fC[0], fC[1]));
ee8fb199 2208 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]));
2209 AliInfo(Form("PID = %5.3f / %5.3f / %5.3f / %5.3f / %5.3f", fProb[0], fProb[1], fProb[2], fProb[3], fProb[4]));
203967fc 2210
2211 if(strcmp(o, "a")!=0) return;
2212
4dc4dc2e 2213 AliTRDcluster* const* jc = &fClusters[0];
8d2bec9e 2214 for(int ic=0; ic<kNclusters; ic++, jc++) {
4dc4dc2e 2215 if(!(*jc)) continue;
203967fc 2216 (*jc)->Print(o);
4dc4dc2e 2217 }
e4f2f73d 2218}
47d5d320 2219
203967fc 2220
2221//___________________________________________________________________
2222Bool_t AliTRDseedV1::IsEqual(const TObject *o) const
2223{
2224 // Checks if current instance of the class has the same essential members
2225 // as the given one
2226
2227 if(!o) return kFALSE;
2228 const AliTRDseedV1 *inTracklet = dynamic_cast<const AliTRDseedV1*>(o);
2229 if(!inTracklet) return kFALSE;
2230
2231 for (Int_t i = 0; i < 2; i++){
e3cf3d02 2232 if ( fYref[i] != inTracklet->fYref[i] ) return kFALSE;
2233 if ( fZref[i] != inTracklet->fZref[i] ) return kFALSE;
203967fc 2234 }
2235
4ecadb52 2236 if ( TMath::Abs(fS2Y - inTracklet->fS2Y)>1.e-10 ) return kFALSE;
2237 if ( TMath::Abs(GetTilt() - inTracklet->GetTilt())>1.e-10 ) return kFALSE;
2238 if ( TMath::Abs(GetPadLength() - inTracklet->GetPadLength())>1.e-10 ) return kFALSE;
203967fc 2239
8d2bec9e 2240 for (Int_t i = 0; i < kNclusters; i++){
e3cf3d02 2241// if ( fX[i] != inTracklet->GetX(i) ) return kFALSE;
2242// if ( fY[i] != inTracklet->GetY(i) ) return kFALSE;
2243// if ( fZ[i] != inTracklet->GetZ(i) ) return kFALSE;
2244 if ( fIndexes[i] != inTracklet->fIndexes[i] ) return kFALSE;
203967fc 2245 }
f29f13a6 2246// if ( fUsable != inTracklet->fUsable ) return kFALSE;
203967fc 2247
2248 for (Int_t i=0; i < 2; i++){
e3cf3d02 2249 if ( fYfit[i] != inTracklet->fYfit[i] ) return kFALSE;
2250 if ( fZfit[i] != inTracklet->fZfit[i] ) return kFALSE;
2251 if ( fLabels[i] != inTracklet->fLabels[i] ) return kFALSE;
203967fc 2252 }
2253
e3cf3d02 2254/* if ( fMeanz != inTracklet->GetMeanz() ) return kFALSE;
2255 if ( fZProb != inTracklet->GetZProb() ) return kFALSE;*/
3e778975 2256 if ( fN != inTracklet->fN ) return kFALSE;
2257 //if ( fNUsed != inTracklet->fNUsed ) return kFALSE;
e3cf3d02 2258 //if ( fFreq != inTracklet->GetFreq() ) return kFALSE;
2259 //if ( fNChange != inTracklet->GetNChange() ) return kFALSE;
203967fc 2260
4ecadb52 2261 if ( TMath::Abs(fC[0] - inTracklet->fC[0])>1.e-10 ) return kFALSE;
e3cf3d02 2262 //if ( fCC != inTracklet->GetCC() ) return kFALSE;
4ecadb52 2263 if ( TMath::Abs(fChi2 - inTracklet->fChi2)>1.e-10 ) return kFALSE;
203967fc 2264 // if ( fChi2Z != inTracklet->GetChi2Z() ) return kFALSE;
2265
e3cf3d02 2266 if ( fDet != inTracklet->fDet ) return kFALSE;
4ecadb52 2267 if ( TMath::Abs(fPt - inTracklet->fPt)>1.e-10 ) return kFALSE;
2268 if ( TMath::Abs(fdX - inTracklet->fdX)>1.e-10 ) return kFALSE;
203967fc 2269
8d2bec9e 2270 for (Int_t iCluster = 0; iCluster < kNclusters; iCluster++){
203967fc 2271 AliTRDcluster *curCluster = fClusters[iCluster];
e3cf3d02 2272 AliTRDcluster *inCluster = inTracklet->fClusters[iCluster];
203967fc 2273 if (curCluster && inCluster){
2274 if (! curCluster->IsEqual(inCluster) ) {
2275 curCluster->Print();
2276 inCluster->Print();
2277 return kFALSE;
2278 }
2279 } else {
2280 // if one cluster exists, and corresponding
2281 // in other tracklet doesn't - return kFALSE
2282 if(curCluster || inCluster) return kFALSE;
2283 }
2284 }
2285 return kTRUE;
2286}
5d401b45 2287