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