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