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