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