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