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