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