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