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