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