- removed unused array
[u/mrichter/AliRoot.git] / TRD / AliTRDseedV1.cxx
CommitLineData
e4f2f73d 1/**************************************************************************
29b87567 2* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3* *
4* Author: The ALICE Off-line Project. *
5* Contributors are mentioned in the code where appropriate. *
6* *
7* Permission to use, copy, modify and distribute this software and its *
8* documentation strictly for non-commercial purposes is hereby granted *
9* without fee, provided that the above copyright notice appears in all *
10* copies and that both the copyright notice and this permission notice *
11* appear in the supporting documentation. The authors make no claims *
12* about the suitability of this software for any purpose. It is *
13* provided "as is" without express or implied warranty. *
14**************************************************************************/
e4f2f73d 15
16/* $Id$ */
17
18////////////////////////////////////////////////////////////////////////////
dd8059a8 19////
20// The TRD offline tracklet
21//
22// The running horse of the TRD reconstruction. The following tasks are preformed:
23// 1. Clusters attachment to tracks based on prior information stored at tracklet level (see AttachClusters)
24// 2. Clusters position recalculation based on track information (see GetClusterXY and Fit)
25// 3. Cluster error parametrization recalculation (see Fit)
26// 4. Linear track approximation (Fit)
27// 5. Optimal position (including z estimate for pad row cross tracklets) and covariance matrix of the track fit inside one TRD chamber (Fit)
28// 6. Tilt pad correction and systematic effects (GetCovAt)
29// 7. dEdx calculation (CookdEdx)
30// 8. PID probabilities estimation (CookPID)
31//
e4f2f73d 32// Authors: //
33// Alex Bercuci <A.Bercuci@gsi.de> //
34// Markus Fasel <M.Fasel@gsi.de> //
35// //
36////////////////////////////////////////////////////////////////////////////
37
38#include "TMath.h"
eb38ed55 39#include <TTreeStream.h>
e4f2f73d 40
41#include "AliLog.h"
42#include "AliMathBase.h"
d937ad7a 43#include "AliCDBManager.h"
44#include "AliTracker.h"
e4f2f73d 45
03cef9b2 46#include "AliTRDpadPlane.h"
e4f2f73d 47#include "AliTRDcluster.h"
f3d3af1b 48#include "AliTRDseedV1.h"
49#include "AliTRDtrackV1.h"
e4f2f73d 50#include "AliTRDcalibDB.h"
eb38ed55 51#include "AliTRDchamberTimeBin.h"
52#include "AliTRDtrackingChamber.h"
53#include "AliTRDtrackerV1.h"
e4f2f73d 54#include "AliTRDrecoParam.h"
a076fc2f 55#include "AliTRDCommonParam.h"
d937ad7a 56
0906e73e 57#include "Cal/AliTRDCalPID.h"
d937ad7a 58#include "Cal/AliTRDCalROC.h"
59#include "Cal/AliTRDCalDet.h"
e4f2f73d 60
e4f2f73d 61ClassImp(AliTRDseedV1)
62
63//____________________________________________________________________
ae4e8b84 64AliTRDseedV1::AliTRDseedV1(Int_t det)
3e778975 65 :AliTRDtrackletBase()
4d6aee34 66 ,fkReconstructor(NULL)
67 ,fClusterIter(NULL)
e3cf3d02 68 ,fExB(0.)
69 ,fVD(0.)
70 ,fT0(0.)
71 ,fS2PRF(0.)
72 ,fDiffL(0.)
73 ,fDiffT(0.)
ae4e8b84 74 ,fClusterIdx(0)
7c3eecb8 75 ,fErrorMsg(0)
3e778975 76 ,fN(0)
ae4e8b84 77 ,fDet(det)
b25a5e9e 78 ,fPt(0.)
bcb6fb78 79 ,fdX(0.)
e3cf3d02 80 ,fX0(0.)
81 ,fX(0.)
82 ,fY(0.)
83 ,fZ(0.)
84 ,fS2Y(0.)
85 ,fS2Z(0.)
e3cf3d02 86 ,fChi2(0.)
e4f2f73d 87{
88 //
89 // Constructor
90 //
f301a656 91 memset(fIndexes,0xFF,kNclusters*sizeof(fIndexes[0]));
8d2bec9e 92 memset(fClusters, 0, kNclusters*sizeof(AliTRDcluster*));
dd8059a8 93 memset(fPad, 0, 3*sizeof(Float_t));
e3cf3d02 94 fYref[0] = 0.; fYref[1] = 0.;
95 fZref[0] = 0.; fZref[1] = 0.;
96 fYfit[0] = 0.; fYfit[1] = 0.;
97 fZfit[0] = 0.; fZfit[1] = 0.;
8d2bec9e 98 memset(fdEdx, 0, kNslices*sizeof(Float_t));
29b87567 99 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = -1.;
e3cf3d02 100 fLabels[0]=-1; fLabels[1]=-1; // most freq MC labels
101 fLabels[2]=0; // number of different labels for tracklet
16cca13f 102 memset(fRefCov, 0, 7*sizeof(Double_t));
68f9b6bd 103 // stand alone curvature
104 fC[0] = 0.; fC[1] = 0.;
d937ad7a 105 // covariance matrix [diagonal]
106 // default sy = 200um and sz = 2.3 cm
107 fCov[0] = 4.e-4; fCov[1] = 0.; fCov[2] = 5.3;
f29f13a6 108 SetStandAlone(kFALSE);
e4f2f73d 109}
110
111//____________________________________________________________________
0906e73e 112AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
3e778975 113 :AliTRDtrackletBase((AliTRDtrackletBase&)ref)
4d6aee34 114 ,fkReconstructor(NULL)
115 ,fClusterIter(NULL)
e3cf3d02 116 ,fExB(0.)
117 ,fVD(0.)
118 ,fT0(0.)
119 ,fS2PRF(0.)
120 ,fDiffL(0.)
121 ,fDiffT(0.)
ae4e8b84 122 ,fClusterIdx(0)
7c3eecb8 123 ,fErrorMsg(0)
3e778975 124 ,fN(0)
e3cf3d02 125 ,fDet(-1)
b25a5e9e 126 ,fPt(0.)
e3cf3d02 127 ,fdX(0.)
128 ,fX0(0.)
129 ,fX(0.)
130 ,fY(0.)
131 ,fZ(0.)
132 ,fS2Y(0.)
133 ,fS2Z(0.)
e3cf3d02 134 ,fChi2(0.)
e4f2f73d 135{
136 //
137 // Copy Constructor performing a deep copy
138 //
e3cf3d02 139 if(this != &ref){
140 ref.Copy(*this);
141 }
29b87567 142 SetBit(kOwner, kFALSE);
f29f13a6 143 SetStandAlone(ref.IsStandAlone());
fbb2ea06 144}
d9950a5a 145
0906e73e 146
e4f2f73d 147//____________________________________________________________________
148AliTRDseedV1& AliTRDseedV1::operator=(const AliTRDseedV1 &ref)
149{
150 //
151 // Assignment Operator using the copy function
152 //
153
29b87567 154 if(this != &ref){
155 ref.Copy(*this);
156 }
221ab7e0 157 SetBit(kOwner, kFALSE);
158
29b87567 159 return *this;
e4f2f73d 160}
161
162//____________________________________________________________________
163AliTRDseedV1::~AliTRDseedV1()
164{
165 //
166 // Destructor. The RecoParam object belongs to the underlying tracker.
167 //
168
29b87567 169 //printf("I-AliTRDseedV1::~AliTRDseedV1() : Owner[%s]\n", IsOwner()?"YES":"NO");
e4f2f73d 170
e3cf3d02 171 if(IsOwner()) {
8d2bec9e 172 for(int itb=0; itb<kNclusters; itb++){
29b87567 173 if(!fClusters[itb]) continue;
174 //AliInfo(Form("deleting c %p @ %d", fClusters[itb], itb));
175 delete fClusters[itb];
4d6aee34 176 fClusters[itb] = NULL;
29b87567 177 }
e3cf3d02 178 }
e4f2f73d 179}
180
181//____________________________________________________________________
182void AliTRDseedV1::Copy(TObject &ref) const
183{
184 //
185 // Copy function
186 //
187
29b87567 188 //AliInfo("");
189 AliTRDseedV1 &target = (AliTRDseedV1 &)ref;
190
4d6aee34 191 target.fkReconstructor = fkReconstructor;
192 target.fClusterIter = NULL;
e3cf3d02 193 target.fExB = fExB;
194 target.fVD = fVD;
195 target.fT0 = fT0;
196 target.fS2PRF = fS2PRF;
197 target.fDiffL = fDiffL;
198 target.fDiffT = fDiffT;
ae4e8b84 199 target.fClusterIdx = 0;
7c3eecb8 200 target.fErrorMsg = fErrorMsg;
3e778975 201 target.fN = fN;
ae4e8b84 202 target.fDet = fDet;
b25a5e9e 203 target.fPt = fPt;
29b87567 204 target.fdX = fdX;
e3cf3d02 205 target.fX0 = fX0;
206 target.fX = fX;
207 target.fY = fY;
208 target.fZ = fZ;
209 target.fS2Y = fS2Y;
210 target.fS2Z = fS2Z;
e3cf3d02 211 target.fChi2 = fChi2;
29b87567 212
8d2bec9e 213 memcpy(target.fIndexes, fIndexes, kNclusters*sizeof(Int_t));
214 memcpy(target.fClusters, fClusters, kNclusters*sizeof(AliTRDcluster*));
dd8059a8 215 memcpy(target.fPad, fPad, 3*sizeof(Float_t));
e3cf3d02 216 target.fYref[0] = fYref[0]; target.fYref[1] = fYref[1];
217 target.fZref[0] = fZref[0]; target.fZref[1] = fZref[1];
218 target.fYfit[0] = fYfit[0]; target.fYfit[1] = fYfit[1];
219 target.fZfit[0] = fZfit[0]; target.fZfit[1] = fZfit[1];
8d2bec9e 220 memcpy(target.fdEdx, fdEdx, kNslices*sizeof(Float_t));
e3cf3d02 221 memcpy(target.fProb, fProb, AliPID::kSPECIES*sizeof(Float_t));
222 memcpy(target.fLabels, fLabels, 3*sizeof(Int_t));
16cca13f 223 memcpy(target.fRefCov, fRefCov, 7*sizeof(Double_t));
68f9b6bd 224 target.fC[0] = fC[0]; target.fC[1] = fC[1];
e3cf3d02 225 memcpy(target.fCov, fCov, 3*sizeof(Double_t));
29b87567 226
e3cf3d02 227 TObject::Copy(ref);
e4f2f73d 228}
229
0906e73e 230
231//____________________________________________________________
f3d3af1b 232Bool_t AliTRDseedV1::Init(AliTRDtrackV1 *track)
0906e73e 233{
234// Initialize this tracklet using the track information
235//
236// Parameters:
237// track - the TRD track used to initialize the tracklet
238//
239// Detailed description
240// The function sets the starting point and direction of the
241// tracklet according to the information from the TRD track.
242//
243// Caution
244// The TRD track has to be propagated to the beginning of the
245// chamber where the tracklet will be constructed
246//
247
29b87567 248 Double_t y, z;
249 if(!track->GetProlongation(fX0, y, z)) return kFALSE;
16cca13f 250 Update(track);
29b87567 251 return kTRUE;
0906e73e 252}
253
bcb6fb78 254
e3cf3d02 255//_____________________________________________________________________________
980d5a2a 256void AliTRDseedV1::Reset(Option_t *opt)
e3cf3d02 257{
980d5a2a 258//
259// Reset seed. If option opt="c" is given only cluster arrays are cleared.
260//
261 for(Int_t ic=kNclusters; ic--;) fIndexes[ic] = -1;
262 memset(fClusters, 0, kNclusters*sizeof(AliTRDcluster*));
560e5c05 263 fN=0; SetBit(kRowCross, kFALSE);
980d5a2a 264 if(strcmp(opt, "c")==0) return;
265
e3cf3d02 266 fExB=0.;fVD=0.;fT0=0.;fS2PRF=0.;
267 fDiffL=0.;fDiffT=0.;
3e778975 268 fClusterIdx=0;
7c3eecb8 269 fErrorMsg = 0;
dd8059a8 270 fDet=-1;
b25a5e9e 271 fPt=0.;
e3cf3d02 272 fdX=0.;fX0=0.; fX=0.; fY=0.; fZ=0.;
273 fS2Y=0.; fS2Z=0.;
68f9b6bd 274 fC[0]=0.; fC[1]=0.;
275 fChi2 = 0.;
e3cf3d02 276
dd8059a8 277 memset(fPad, 0, 3*sizeof(Float_t));
e3cf3d02 278 fYref[0] = 0.; fYref[1] = 0.;
279 fZref[0] = 0.; fZref[1] = 0.;
280 fYfit[0] = 0.; fYfit[1] = 0.;
281 fZfit[0] = 0.; fZfit[1] = 0.;
8d2bec9e 282 memset(fdEdx, 0, kNslices*sizeof(Float_t));
e3cf3d02 283 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = -1.;
284 fLabels[0]=-1; fLabels[1]=-1; // most freq MC labels
285 fLabels[2]=0; // number of different labels for tracklet
16cca13f 286 memset(fRefCov, 0, 7*sizeof(Double_t));
e3cf3d02 287 // covariance matrix [diagonal]
288 // default sy = 200um and sz = 2.3 cm
289 fCov[0] = 4.e-4; fCov[1] = 0.; fCov[2] = 5.3;
290}
291
bcb6fb78 292//____________________________________________________________________
16cca13f 293void AliTRDseedV1::Update(const AliTRDtrackV1 *trk)
b1957d3c 294{
295 // update tracklet reference position from the TRD track
b1957d3c 296
e3cf3d02 297 Double_t fSnp = trk->GetSnp();
298 Double_t fTgl = trk->GetTgl();
b25a5e9e 299 fPt = trk->Pt();
bfd20868 300 Double_t norm =1./TMath::Sqrt((1.-fSnp)*(1.+fSnp));
1fd9389f 301 fYref[1] = fSnp*norm;
302 fZref[1] = fTgl*norm;
b1957d3c 303 SetCovRef(trk->GetCovariance());
304
305 Double_t dx = trk->GetX() - fX0;
306 fYref[0] = trk->GetY() - dx*fYref[1];
307 fZref[0] = trk->GetZ() - dx*fZref[1];
308}
309
e3cf3d02 310//_____________________________________________________________________________
311void AliTRDseedV1::UpdateUsed()
312{
313 //
f29f13a6 314 // Calculate number of used clusers in the tracklet
e3cf3d02 315 //
316
3e778975 317 Int_t nused = 0, nshared = 0;
8d2bec9e 318 for (Int_t i = kNclusters; i--; ) {
e3cf3d02 319 if (!fClusters[i]) continue;
3e778975 320 if(fClusters[i]->IsUsed()){
321 nused++;
322 } else if(fClusters[i]->IsShared()){
323 if(IsStandAlone()) nused++;
324 else nshared++;
325 }
e3cf3d02 326 }
3e778975 327 SetNUsed(nused);
328 SetNShared(nshared);
e3cf3d02 329}
330
331//_____________________________________________________________________________
332void AliTRDseedV1::UseClusters()
333{
334 //
335 // Use clusters
336 //
f29f13a6 337 // In stand alone mode:
338 // Clusters which are marked as used or shared from another track are
339 // removed from the tracklet
340 //
341 // In barrel mode:
342 // - Clusters which are used by another track become shared
343 // - Clusters which are attached to a kink track become shared
344 //
e3cf3d02 345 AliTRDcluster **c = &fClusters[0];
8d2bec9e 346 for (Int_t ic=kNclusters; ic--; c++) {
e3cf3d02 347 if(!(*c)) continue;
f29f13a6 348 if(IsStandAlone()){
349 if((*c)->IsShared() || (*c)->IsUsed()){
b82b4de1 350 if((*c)->IsShared()) SetNShared(GetNShared()-1);
351 else SetNUsed(GetNUsed()-1);
4d6aee34 352 (*c) = NULL;
f29f13a6 353 fIndexes[ic] = -1;
3e778975 354 SetN(GetN()-1);
3e778975 355 continue;
f29f13a6 356 }
3e778975 357 } else {
f29f13a6 358 if((*c)->IsUsed() || IsKink()){
3e778975 359 (*c)->SetShared();
360 continue;
f29f13a6 361 }
362 }
363 (*c)->Use();
e3cf3d02 364 }
365}
366
367
f29f13a6 368
b1957d3c 369//____________________________________________________________________
bcb6fb78 370void AliTRDseedV1::CookdEdx(Int_t nslices)
371{
372// Calculates average dE/dx for all slices and store them in the internal array fdEdx.
373//
374// Parameters:
375// nslices : number of slices for which dE/dx should be calculated
376// Output:
377// store results in the internal array fdEdx. This can be accessed with the method
378// AliTRDseedV1::GetdEdx()
379//
380// Detailed description
381// Calculates average dE/dx for all slices. Depending on the PID methode
382// the number of slices can be 3 (LQ) or 8(NN).
3ee48d6e 383// The calculation of dQ/dl are done using the tracklet fit results (see AliTRDseedV1::GetdQdl(Int_t))
bcb6fb78 384//
385// The following effects are included in the calculation:
386// 1. calibration values for t0 and vdrift (using x coordinate to calculate slice)
387// 2. cluster sharing (optional see AliTRDrecoParam::SetClusterSharing())
388// 3. cluster size
389//
390
8d2bec9e 391 memset(fdEdx, 0, kNslices*sizeof(Float_t));
e73abf77 392 const Double_t kDriftLength = (.5 * AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
29b87567 393
9ded305e 394 AliTRDcluster *c(NULL);
29b87567 395 for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++){
8e709c82 396 if(!(c = fClusters[ic]) && !(c = fClusters[ic+kNtb])) continue;
e73abf77 397 Float_t dx = TMath::Abs(fX0 - c->GetX());
29b87567 398
399 // Filter clusters for dE/dx calculation
400
401 // 1.consider calibration effects for slice determination
e73abf77 402 Int_t slice;
403 if(dx<kDriftLength){ // TODO should be replaced by c->IsInChamber()
404 slice = Int_t(dx * nslices / kDriftLength);
405 } else slice = c->GetX() < fX0 ? nslices-1 : 0;
406
407
29b87567 408 // 2. take sharing into account
3e778975 409 Float_t w = /*c->IsShared() ? .5 :*/ 1.;
29b87567 410
411 // 3. take into account large clusters TODO
412 //w *= c->GetNPads() > 3 ? .8 : 1.;
413
414 //CHECK !!!
415 fdEdx[slice] += w * GetdQdl(ic); //fdQdl[ic];
29b87567 416 } // End of loop over clusters
bcb6fb78 417}
418
e3cf3d02 419//_____________________________________________________________________________
420void AliTRDseedV1::CookLabels()
421{
422 //
423 // Cook 2 labels for seed
424 //
425
426 Int_t labels[200];
427 Int_t out[200];
428 Int_t nlab = 0;
8d2bec9e 429 for (Int_t i = 0; i < kNclusters; i++) {
e3cf3d02 430 if (!fClusters[i]) continue;
431 for (Int_t ilab = 0; ilab < 3; ilab++) {
432 if (fClusters[i]->GetLabel(ilab) >= 0) {
433 labels[nlab] = fClusters[i]->GetLabel(ilab);
434 nlab++;
435 }
436 }
437 }
438
fac58f00 439 fLabels[2] = AliMathBase::Freq(nlab,labels,out,kTRUE);
e3cf3d02 440 fLabels[0] = out[0];
441 if ((fLabels[2] > 1) && (out[3] > 1)) fLabels[1] = out[2];
442}
443
444
d937ad7a 445//____________________________________________________________________
0b433f72 446Float_t AliTRDseedV1::GetdQdl(Int_t ic, Float_t *dl) const
bcb6fb78 447{
3ee48d6e 448// Using the linear approximation of the track inside one TRD chamber (TRD tracklet)
449// the charge per unit length can be written as:
450// BEGIN_LATEX
500851ab 451// #frac{dq}{dl} = #frac{q_{c}}{dx * #sqrt{1 + #(){#frac{dy}{dx}}^{2}_{fit} + #(){#frac{dz}{dx}}^{2}_{ref}}}
3ee48d6e 452// END_LATEX
453// where qc is the total charge collected in the current time bin and dx is the length
0b433f72 454// of the time bin.
455// The following correction are applied :
456// - charge : pad row cross corrections
457// [diffusion and TRF assymetry] TODO
458// - dx : anisochronity, track inclination - see Fit and AliTRDcluster::GetXloc()
459// and AliTRDcluster::GetYloc() for the effects taken into account
3ee48d6e 460//
0fa1a8ee 461//Begin_Html
462//<img src="TRD/trackletDQDT.gif">
463//End_Html
464// In the picture the energy loss measured on the tracklet as a function of drift time [left] and respectively
465// drift length [right] for different particle species is displayed.
3ee48d6e 466// Author : Alex Bercuci <A.Bercuci@gsi.de>
467//
468 Float_t dq = 0.;
5d401b45 469 // check whether both clusters are inside the chamber
470 Bool_t hasClusterInChamber = kFALSE;
471 if(fClusters[ic] && fClusters[ic]->IsInChamber()){
472 hasClusterInChamber = kTRUE;
1742f24c 473 dq += TMath::Abs(fClusters[ic]->GetQ());
5d401b45 474 }else if(fClusters[ic+kNtb] && fClusters[ic+kNtb]->IsInChamber()){
475 hasClusterInChamber = kTRUE;
476 dq += TMath::Abs(fClusters[ic+kNtb]->GetQ());
1742f24c 477 }
5d401b45 478 if(!hasClusterInChamber) return 0.;
0b433f72 479 if(dq<1.e-3) return 0.;
3ee48d6e 480
a2abcbc5 481 Double_t dx = fdX;
482 if(ic-1>=0 && ic+1<kNtb){
483 Float_t x2(0.), x1(0.);
5d401b45 484 // try to estimate upper radial position (find the cluster which is inside the chamber)
485 if(fClusters[ic-1] && fClusters[ic-1]->IsInChamber()) x2 = fClusters[ic-1]->GetX();
486 else if(fClusters[ic-1+kNtb] && fClusters[ic-1+kNtb]->IsInChamber()) x2 = fClusters[ic-1+kNtb]->GetX();
487 else if(fClusters[ic] && fClusters[ic]->IsInChamber()) x2 = fClusters[ic]->GetX()+fdX;
a2abcbc5 488 else x2 = fClusters[ic+kNtb]->GetX()+fdX;
5d401b45 489 // try to estimate lower radial position (find the cluster which is inside the chamber)
490 if(fClusters[ic+1] && fClusters[ic+1]->IsInChamber()) x1 = fClusters[ic+1]->GetX();
491 else if(fClusters[ic+1+kNtb] && fClusters[ic+1+kNtb]->IsInChamber()) x1 = fClusters[ic+1+kNtb]->GetX();
492 else if(fClusters[ic] && fClusters[ic]->IsInChamber()) x1 = fClusters[ic]->GetX()-fdX;
a2abcbc5 493 else x1 = fClusters[ic+kNtb]->GetX()-fdX;
494
495 dx = .5*(x2 - x1);
496 }
0b433f72 497 dx *= TMath::Sqrt(1. + fYfit[1]*fYfit[1] + fZref[1]*fZref[1]);
0b433f72 498 if(dl) (*dl) = dx;
283604d2 499 if(dx>1.e-9) return dq/dx;
500 else return 0.;
bcb6fb78 501}
502
0b433f72 503//____________________________________________________________
504Float_t AliTRDseedV1::GetMomentum(Float_t *err) const
505{
506// Returns momentum of the track after update with the current tracklet as:
507// BEGIN_LATEX
508// p=#frac{1}{1/p_{t}} #sqrt{1+tgl^{2}}
509// END_LATEX
510// and optionally the momentum error (if err is not null).
511// The estimated variance of the momentum is given by:
512// BEGIN_LATEX
513// #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})
514// END_LATEX
515// which can be simplified to
516// BEGIN_LATEX
517// #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}
518// END_LATEX
519//
520
521 Double_t p = fPt*TMath::Sqrt(1.+fZref[1]*fZref[1]);
522 Double_t p2 = p*p;
523 Double_t tgl2 = fZref[1]*fZref[1];
524 Double_t pt2 = fPt*fPt;
525 if(err){
526 Double_t s2 =
527 p2*tgl2*pt2*pt2*fRefCov[4]
528 -2.*p2*fZref[1]*fPt*pt2*fRefCov[5]
529 +p2*pt2*fRefCov[6];
530 (*err) = TMath::Sqrt(s2);
531 }
532 return p;
533}
534
b453ef55 535//____________________________________________________________________
536Float_t AliTRDseedV1::GetOccupancyTB() const
537{
538// Returns procentage of TB occupied by clusters
539
540 Int_t n(0);
541 AliTRDcluster *c(NULL);
542 for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++){
543 if(!(c = fClusters[ic]) && !(c = fClusters[ic+kNtb])) continue;
544 n++;
545 }
546
547 return Float_t(n)/AliTRDtrackerV1::GetNTimeBins();
548}
0b433f72 549
0906e73e 550//____________________________________________________________________
3e778975 551Float_t* AliTRDseedV1::GetProbability(Bool_t force)
0906e73e 552{
3e778975 553 if(!force) return &fProb[0];
4d6aee34 554 if(!CookPID()) return NULL;
3e778975 555 return &fProb[0];
556}
557
558//____________________________________________________________
559Bool_t AliTRDseedV1::CookPID()
560{
0906e73e 561// Fill probability array for tracklet from the DB.
562//
563// Parameters
564//
565// Output
4d6aee34 566// returns pointer to the probability array and NULL if missing DB access
0906e73e 567//
2a3191bb 568// Retrieve PID probabilities for e+-, mu+-, K+-, pi+- and p+- from the DB according to tracklet information:
569// - estimated momentum at tracklet reference point
570// - dE/dx measurements
571// - tracklet length
572// - TRD layer
573// According to the steering settings specified in the reconstruction one of the following methods are used
574// - Neural Network [default] - option "nn"
575// - 2D Likelihood - option "!nn"
0906e73e 576
0906e73e 577 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
578 if (!calibration) {
579 AliError("No access to calibration data");
3e778975 580 return kFALSE;
0906e73e 581 }
582
4d6aee34 583 if (!fkReconstructor) {
3a039a31 584 AliError("Reconstructor not set.");
3e778975 585 return kFALSE;
4ba1d6ae 586 }
587
0906e73e 588 // Retrieve the CDB container class with the parametric detector response
4d6aee34 589 const AliTRDCalPID *pd = calibration->GetPIDObject(fkReconstructor->GetPIDMethod());
0906e73e 590 if (!pd) {
591 AliError("No access to AliTRDCalPID object");
3e778975 592 return kFALSE;
0906e73e 593 }
10f75631 594
29b87567 595 // calculate tracklet length TO DO
560e5c05 596 Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())/ TMath::Sqrt((1.0 - GetSnp()*GetSnp()) / (1.0 + GetTgl()*GetTgl()));
0906e73e 597
598 //calculate dE/dx
9ded305e 599 CookdEdx(AliTRDCalPID::kNSlicesNN);
600 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 601
0906e73e 602 // Sets the a priori probabilities
11d80e40 603 Bool_t kPIDNN(fkReconstructor->GetPIDMethod()==AliTRDpidUtil::kNN);
f83cd814 604 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++)
11d80e40 605 fProb[ispec] = pd->GetProbability(ispec, GetMomentum(), &fdEdx[0], length, kPIDNN?GetPlane():fkReconstructor->GetRecoParam()->GetPIDLQslices());
f301a656 606
3e778975 607 return kTRUE;
0906e73e 608}
609
610//____________________________________________________________________
e4f2f73d 611Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
612{
613 //
614 // Returns a quality measurement of the current seed
615 //
616
dd8059a8 617 Float_t zcorr = kZcorr ? GetTilt() * (fZfit[0] - fZref[0]) : 0.;
29b87567 618 return
3e778975 619 .5 * TMath::Abs(18.0 - GetN())
29b87567 620 + 10.* TMath::Abs(fYfit[1] - fYref[1])
621 + 5. * TMath::Abs(fYfit[0] - fYref[0] + zcorr)
dd8059a8 622 + 2. * TMath::Abs(fZfit[0] - fZref[0]) / GetPadLength();
e4f2f73d 623}
624
625//____________________________________________________________________
d937ad7a 626void AliTRDseedV1::GetCovAt(Double_t x, Double_t *cov) const
0906e73e 627{
d937ad7a 628// Computes covariance in the y-z plane at radial point x (in tracking coordinates)
629// and returns the results in the preallocated array cov[3] as :
630// cov[0] = Var(y)
631// cov[1] = Cov(yz)
632// cov[2] = Var(z)
633//
634// Details
635//
636// For the linear transformation
637// BEGIN_LATEX
638// Y = T_{x} X^{T}
639// END_LATEX
640// The error propagation has the general form
641// BEGIN_LATEX
642// C_{Y} = T_{x} C_{X} T_{x}^{T}
643// END_LATEX
644// We apply this formula 2 times. First to calculate the covariance of the tracklet
645// at point x we consider:
646// BEGIN_LATEX
647// T_{x} = (1 x); X=(y0 dy/dx); C_{X}=#(){#splitline{Var(y0) Cov(y0, dy/dx)}{Cov(y0, dy/dx) Var(dy/dx)}}
648// END_LATEX
649// and secondly to take into account the tilt angle
650// BEGIN_LATEX
651// T_{#alpha} = #(){#splitline{cos(#alpha) __ sin(#alpha)}{-sin(#alpha) __ cos(#alpha)}}; X=(y z); C_{X}=#(){#splitline{Var(y) 0}{0 Var(z)}}
652// END_LATEX
653//
654// using simple trigonometrics one can write for this last case
655// BEGIN_LATEX
656// 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})}}
657// END_LATEX
658// which can be aproximated for small alphas (2 deg) with
659// BEGIN_LATEX
660// 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}}}
661// END_LATEX
662//
663// before applying the tilt rotation we also apply systematic uncertainties to the tracklet
664// position which can be tunned from outside via the AliTRDrecoParam::SetSysCovMatrix(). They might
665// account for extra misalignment/miscalibration uncertainties.
666//
667// Author :
668// Alex Bercuci <A.Bercuci@gsi.de>
669// Date : Jan 8th 2009
670//
b1957d3c 671
672
d937ad7a 673 Double_t xr = fX0-x;
674 Double_t sy2 = fCov[0] +2.*xr*fCov[1] + xr*xr*fCov[2];
b72f4eaf 675 Double_t sz2 = fS2Z;
676 //GetPadLength()*GetPadLength()/12.;
0906e73e 677
d937ad7a 678 // insert systematic uncertainties
4d6aee34 679 if(fkReconstructor){
bb2db46c 680 Double_t sys[15]; memset(sys, 0, 15*sizeof(Double_t));
4d6aee34 681 fkReconstructor->GetRecoParam()->GetSysCovMatrix(sys);
bb2db46c 682 sy2 += sys[0];
683 sz2 += sys[1];
684 }
d937ad7a 685 // rotate covariance matrix
dd8059a8 686 Double_t t2 = GetTilt()*GetTilt();
d937ad7a 687 Double_t correction = 1./(1. + t2);
688 cov[0] = (sy2+t2*sz2)*correction;
dd8059a8 689 cov[1] = GetTilt()*(sz2 - sy2)*correction;
d937ad7a 690 cov[2] = (t2*sy2+sz2)*correction;
b72f4eaf 691
692 //printf("C(%6.1f %+6.3f %6.1f) [%s]\n", 1.e4*TMath::Sqrt(cov[0]), cov[1], 1.e4*TMath::Sqrt(cov[2]), IsRowCross()?" RC ":"-");
d937ad7a 693}
eb38ed55 694
bb2db46c 695//____________________________________________________________
66765e8e 696Int_t AliTRDseedV1::GetCovSqrt(const Double_t * const c, Double_t *d)
bb2db46c 697{
698// Helper function to calculate the square root of the covariance matrix.
699// The input matrix is stored in the vector c and the result in the vector d.
41b7c7b6 700// Both arrays have to be initialized by the user with at least 3 elements. Return negative in case of failure.
bb2db46c 701//
ec3f0161 702// For calculating the square root of the symmetric matrix c
703// the following relation is used:
bb2db46c 704// BEGIN_LATEX
ec3f0161 705// C^{1/2} = VD^{1/2}V^{-1}
bb2db46c 706// END_LATEX
41b7c7b6 707// with V being the matrix with the n eigenvectors as columns.
ec3f0161 708// In case C is symmetric the followings are true:
709// - matrix D is diagonal with the diagonal given by the eigenvalues of C
41b7c7b6 710// - V = V^{-1}
bb2db46c 711//
712// Author A.Bercuci <A.Bercuci@gsi.de>
713// Date Mar 19 2009
714
66765e8e 715 const Double_t kZero(1.e-20);
4d6aee34 716 Double_t l[2], // eigenvalues
717 v[3]; // eigenvectors
bb2db46c 718 // the secular equation and its solution :
719 // (c[0]-L)(c[2]-L)-c[1]^2 = 0
720 // L^2 - L*Tr(c)+DET(c) = 0
721 // L12 = [Tr(c) +- sqrt(Tr(c)^2-4*DET(c))]/2
4d6aee34 722 Double_t tr = c[0]+c[2], // trace
723 det = c[0]*c[2]-c[1]*c[1]; // determinant
66765e8e 724 if(TMath::Abs(det)<kZero) return 1;
4d6aee34 725 Double_t dd = TMath::Sqrt(tr*tr - 4*det);
66765e8e 726 l[0] = .5*(tr + dd*(c[0]>c[2]?-1.:1.));
727 l[1] = .5*(tr + dd*(c[0]>c[2]?1.:-1.));
728 if(l[0]<kZero || l[1]<kZero) return 2;
41b7c7b6 729 // the sym V matrix
730 // | v00 v10|
731 // | v10 v11|
66765e8e 732 Double_t den = (l[0]-c[0])*(l[0]-c[0])+c[1]*c[1];
733 if(den<kZero){ // almost diagonal
734 v[0] = TMath::Sign(0., c[1]);
735 v[1] = TMath::Sign(1., (l[0]-c[0]));
736 v[2] = TMath::Sign(0., c[1]*(l[0]-c[0])*(l[1]-c[2]));
737 } else {
738 Double_t tmp = 1./TMath::Sqrt(den);
739 v[0] = c[1]* tmp;
740 v[1] = (l[0]-c[0])*tmp;
741 if(TMath::Abs(l[1]-c[2])<kZero) v[2] = TMath::Sign(v[0]*(l[0]-c[0])/kZero, (l[1]-c[2]));
742 else v[2] = v[0]*(l[0]-c[0])/(l[1]-c[2]);
743 }
41b7c7b6 744 // the VD^{1/2}V is:
4d6aee34 745 l[0] = TMath::Sqrt(l[0]); l[1] = TMath::Sqrt(l[1]);
746 d[0] = v[0]*v[0]*l[0]+v[1]*v[1]*l[1];
747 d[1] = v[0]*v[1]*l[0]+v[1]*v[2]*l[1];
748 d[2] = v[1]*v[1]*l[0]+v[2]*v[2]*l[1];
bb2db46c 749
66765e8e 750 return 0;
bb2db46c 751}
752
753//____________________________________________________________
4d6aee34 754Double_t AliTRDseedV1::GetCovInv(const Double_t * const c, Double_t *d)
bb2db46c 755{
756// Helper function to calculate the inverse of the covariance matrix.
757// The input matrix is stored in the vector c and the result in the vector d.
758// Both arrays have to be initialized by the user with at least 3 elements
759// The return value is the determinant or 0 in case of singularity.
760//
761// Author A.Bercuci <A.Bercuci@gsi.de>
762// Date Mar 19 2009
763
4d6aee34 764 Double_t det = c[0]*c[2] - c[1]*c[1];
765 if(TMath::Abs(det)<1.e-20) return 0.;
766 Double_t invDet = 1./det;
767 d[0] = c[2]*invDet;
768 d[1] =-c[1]*invDet;
769 d[2] = c[0]*invDet;
770 return det;
bb2db46c 771}
0906e73e 772
d937ad7a 773//____________________________________________________________________
b72f4eaf 774UShort_t AliTRDseedV1::GetVolumeId() const
775{
fbe11be7 776 for(Int_t ic(0);ic<kNclusters; ic++){
777 if(fClusters[ic]) return fClusters[ic]->GetVolumeId();
778 }
779 return 0;
b72f4eaf 780}
781
782
783//____________________________________________________________________
e3cf3d02 784void AliTRDseedV1::Calibrate()
d937ad7a 785{
e3cf3d02 786// Retrieve calibration and position parameters from OCDB.
787// The following information are used
d937ad7a 788// - detector index
e3cf3d02 789// - column and row position of first attached cluster. If no clusters are attached
790// to the tracklet a random central chamber position (c=70, r=7) will be used.
791//
792// The following information is cached in the tracklet
793// t0 (trigger delay)
794// drift velocity
795// PRF width
796// omega*tau = tg(a_L)
797// diffusion coefficients (longitudinal and transversal)
d937ad7a 798//
799// Author :
800// Alex Bercuci <A.Bercuci@gsi.de>
801// Date : Jan 8th 2009
802//
eb38ed55 803
d937ad7a 804 AliCDBManager *cdb = AliCDBManager::Instance();
805 if(cdb->GetRun() < 0){
806 AliError("OCDB manager not properly initialized");
807 return;
808 }
0906e73e 809
e3cf3d02 810 AliTRDcalibDB *calib = AliTRDcalibDB::Instance();
811 AliTRDCalROC *vdROC = calib->GetVdriftROC(fDet),
812 *t0ROC = calib->GetT0ROC(fDet);;
813 const AliTRDCalDet *vdDet = calib->GetVdriftDet();
814 const AliTRDCalDet *t0Det = calib->GetT0Det();
d937ad7a 815
816 Int_t col = 70, row = 7;
817 AliTRDcluster **c = &fClusters[0];
3e778975 818 if(GetN()){
d937ad7a 819 Int_t ic = 0;
8d2bec9e 820 while (ic<kNclusters && !(*c)){ic++; c++;}
d937ad7a 821 if(*c){
822 col = (*c)->GetPadCol();
823 row = (*c)->GetPadRow();
824 }
825 }
3a039a31 826
e17f4785 827 fT0 = (t0Det->GetValue(fDet) + t0ROC->GetValue(col,row)) / AliTRDCommonParam::Instance()->GetSamplingFrequency();
e3cf3d02 828 fVD = vdDet->GetValue(fDet) * vdROC->GetValue(col, row);
829 fS2PRF = calib->GetPRFWidth(fDet, col, row); fS2PRF *= fS2PRF;
830 fExB = AliTRDCommonParam::Instance()->GetOmegaTau(fVD);
831 AliTRDCommonParam::Instance()->GetDiffCoeff(fDiffL,
832 fDiffT, fVD);
903326c1 833 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));
834
835
e3cf3d02 836 SetBit(kCalib, kTRUE);
0906e73e 837}
838
0906e73e 839//____________________________________________________________________
29b87567 840void AliTRDseedV1::SetOwner()
0906e73e 841{
29b87567 842 //AliInfo(Form("own [%s] fOwner[%s]", own?"YES":"NO", fOwner?"YES":"NO"));
843
844 if(TestBit(kOwner)) return;
8d2bec9e 845 for(int ic=0; ic<kNclusters; ic++){
29b87567 846 if(!fClusters[ic]) continue;
847 fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
848 }
849 SetBit(kOwner);
0906e73e 850}
851
eb2b4f91 852//____________________________________________________________
853void AliTRDseedV1::SetPadPlane(AliTRDpadPlane *p)
854{
855// Shortcut method to initialize pad geometry.
856 if(!p) return;
857 SetTilt(TMath::Tan(TMath::DegToRad()*p->GetTiltingAngle()));
858 SetPadLength(p->GetLengthIPad());
859 SetPadWidth(p->GetWidthIPad());
860}
861
862
e4f2f73d 863//____________________________________________________________________
4d6aee34 864Bool_t AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *const chamber, Bool_t tilt)
e4f2f73d 865{
1fd9389f 866//
867// Projective algorithm to attach clusters to seeding tracklets. The following steps are performed :
868// 1. Collapse x coordinate for the full detector plane
869// 2. truncated mean on y (r-phi) direction
870// 3. purge clusters
871// 4. truncated mean on z direction
872// 5. purge clusters
873//
874// Parameters
875// - chamber : pointer to tracking chamber container used to search the tracklet
876// - tilt : switch for tilt correction during road building [default true]
877// Output
878// - true : if tracklet found successfully. Failure can happend because of the following:
879// -
880// Detailed description
881//
882// We start up by defining the track direction in the xy plane and roads. The roads are calculated based
8a7ff53c 883// on tracking information (variance in the r-phi direction) and estimated variance of the standard
884// clusters (see AliTRDcluster::SetSigmaY2()) corrected for tilt (see GetCovAt()). From this the road is
885// BEGIN_LATEX
500851ab 886// 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 887// r_{z} = 1.5*L_{pad}
888// END_LATEX
1fd9389f 889//
4b755889 890// Author : Alexandru Bercuci <A.Bercuci@gsi.de>
891// Debug : level >3
1fd9389f 892
fc0882f3 893 const AliTRDrecoParam* const recoParam = fkReconstructor->GetRecoParam(); //the dynamic cast in GetRecoParam is slow, so caching the pointer to it
894
895 if(!recoParam){
560e5c05 896 AliError("Tracklets can not be used without a valid RecoParam.");
29b87567 897 return kFALSE;
898 }
b1957d3c 899 // Initialize reco params for this tracklet
900 // 1. first time bin in the drift region
a2abcbc5 901 Int_t t0 = 14;
fc0882f3 902 Int_t kClmin = Int_t(recoParam->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
29b87567 903
fc0882f3 904 Double_t sysCov[5]; recoParam->GetSysCovMatrix(sysCov);
8a7ff53c 905 Double_t s2yTrk= fRefCov[0],
906 s2yCl = 0.,
907 s2zCl = GetPadLength()*GetPadLength()/12.,
908 syRef = TMath::Sqrt(s2yTrk),
909 t2 = GetTilt()*GetTilt();
29b87567 910 //define roads
fc0882f3 911 Double_t kroady = 1., //recoParam->GetRoad1y();
912 kroadz = GetPadLength() * recoParam->GetRoadzMultiplicator() + 1.;
8a7ff53c 913 // define probing cluster (the perfect cluster) and default calibration
914 Short_t sig[] = {0, 0, 10, 30, 10, 0,0};
915 AliTRDcluster cp(fDet, 6, 75, 0, sig, 0);
560e5c05 916 if(fkReconstructor->IsHLT()) cp.SetRPhiMethod(AliTRDcluster::kCOG);
917 if(!IsCalibrated()) Calibrate();
8a7ff53c 918
ee8fb199 919 AliDebug(4, "");
920 AliDebug(4, Form("syKalman[%f] rY[%f] rZ[%f]", syRef, kroady, kroadz));
29b87567 921
922 // working variables
b1957d3c 923 const Int_t kNrows = 16;
4b755889 924 const Int_t kNcls = 3*kNclusters; // buffer size
925 AliTRDcluster *clst[kNrows][kNcls];
3044dfe5 926 Bool_t blst[kNrows][kNcls];
4b755889 927 Double_t cond[4], dx, dy, yt, zt, yres[kNrows][kNcls];
928 Int_t idxs[kNrows][kNcls], ncl[kNrows], ncls = 0;
b1957d3c 929 memset(ncl, 0, kNrows*sizeof(Int_t));
4b755889 930 memset(yres, 0, kNrows*kNcls*sizeof(Double_t));
3044dfe5 931 memset(blst, 0, kNrows*kNcls*sizeof(Bool_t)); //this is 8 times faster to memset than "memset(clst, 0, kNrows*kNcls*sizeof(AliTRDcluster*))"
b1957d3c 932
29b87567 933 // Do cluster projection
4d6aee34 934 AliTRDcluster *c = NULL;
935 AliTRDchamberTimeBin *layer = NULL;
b1957d3c 936 Bool_t kBUFFER = kFALSE;
4b755889 937 for (Int_t it = 0; it < kNtb; it++) {
b1957d3c 938 if(!(layer = chamber->GetTB(it))) continue;
29b87567 939 if(!Int_t(*layer)) continue;
8a7ff53c 940 // get track projection at layers position
b1957d3c 941 dx = fX0 - layer->GetX();
942 yt = fYref[0] - fYref[1] * dx;
943 zt = fZref[0] - fZref[1] * dx;
8a7ff53c 944 // get standard cluster error corrected for tilt
945 cp.SetLocalTimeBin(it);
946 cp.SetSigmaY2(0.02, fDiffT, fExB, dx, -1./*zt*/, fYref[1]);
d956a643 947 s2yCl = (cp.GetSigmaY2() + sysCov[0] + t2*s2zCl)/(1.+t2);
8a7ff53c 948 // get estimated road
949 kroady = 3.*TMath::Sqrt(12.*(s2yTrk + s2yCl));
950
ee8fb199 951 AliDebug(5, Form(" %2d x[%f] yt[%f] zt[%f]", it, dx, yt, zt));
952
953 AliDebug(5, Form(" syTrk[um]=%6.2f syCl[um]=%6.2f syClTlt[um]=%6.2f Ry[mm]=%f", 1.e4*TMath::Sqrt(s2yTrk), 1.e4*TMath::Sqrt(cp.GetSigmaY2()), 1.e4*TMath::Sqrt(s2yCl), 1.e1*kroady));
b1957d3c 954
8a7ff53c 955 // select clusters
b1957d3c 956 cond[0] = yt; cond[2] = kroady;
957 cond[1] = zt; cond[3] = kroadz;
958 Int_t n=0, idx[6];
959 layer->GetClusters(cond, idx, n, 6);
960 for(Int_t ic = n; ic--;){
961 c = (*layer)[idx[ic]];
962 dy = yt - c->GetY();
dd8059a8 963 dy += tilt ? GetTilt() * (c->GetZ() - zt) : 0.;
b1957d3c 964 // select clusters on a 3 sigmaKalman level
965/* if(tilt && TMath::Abs(dy) > 3.*syRef){
966 printf("too large !!!\n");
967 continue;
968 }*/
969 Int_t r = c->GetPadRow();
ee8fb199 970 AliDebug(5, Form(" -> dy[%f] yc[%f] r[%d]", TMath::Abs(dy), c->GetY(), r));
b1957d3c 971 clst[r][ncl[r]] = c;
3044dfe5 972 blst[r][ncl[r]] = kTRUE;
b1957d3c 973 idxs[r][ncl[r]] = idx[ic];
974 yres[r][ncl[r]] = dy;
975 ncl[r]++; ncls++;
976
4b755889 977 if(ncl[r] >= kNcls) {
560e5c05 978 AliWarning(Form("Cluster candidates row[%d] reached buffer limit[%d]. Some may be lost.", r, kNcls));
b1957d3c 979 kBUFFER = kTRUE;
29b87567 980 break;
981 }
982 }
b1957d3c 983 if(kBUFFER) break;
29b87567 984 }
ee8fb199 985 AliDebug(4, Form("Found %d clusters. Processing ...", ncls));
986 if(ncls<kClmin){
560e5c05 987 AliDebug(1, Form("CLUSTERS FOUND %d LESS THAN THRESHOLD %d.", ncls, kClmin));
7c3eecb8 988 SetErrorMsg(kAttachClFound);
ee8fb199 989 return kFALSE;
990 }
991
b1957d3c 992 // analyze each row individualy
560e5c05 993 Bool_t kRowSelection(kFALSE);
994 Double_t mean[]={1.e3, 1.e3, 1.3}, syDis[]={1.e3, 1.e3, 1.3};
995 Int_t nrow[] = {0, 0, 0}, rowId[] = {-1, -1, -1}, nr = 0, lr=-1;
996 TVectorD vdy[3];
997 for(Int_t ir=0; ir<kNrows; ir++){
b1957d3c 998 if(!(ncl[ir])) continue;
560e5c05 999 if(lr>0 && ir-lr != 1){
1000 AliDebug(2, "Rows attached not continuous. Turn on selection.");
1001 kRowSelection=kTRUE;
29b87567 1002 }
560e5c05 1003
ee8fb199 1004 AliDebug(5, Form(" r[%d] n[%d]", ir, ncl[ir]));
b1957d3c 1005 // Evaluate truncated mean on the y direction
560e5c05 1006 if(ncl[ir] < 4) continue;
1007 AliMathBase::EvaluateUni(ncl[ir], yres[ir], mean[nr], syDis[nr], Int_t(ncl[ir]*.8));
4b755889 1008
b1957d3c 1009 // TODO check mean and sigma agains cluster resolution !!
560e5c05 1010 AliDebug(4, Form(" m_%d[%+5.3f (%5.3fs)] s[%f]", nr, mean[nr], TMath::Abs(mean[nr]/syDis[nr]), syDis[nr]));
1011 // remove outliers based on a 3 sigmaDistr level
b1957d3c 1012 Bool_t kFOUND = kFALSE;
1013 for(Int_t ic = ncl[ir]; ic--;){
560e5c05 1014 if(yres[ir][ic] - mean[nr] > 3. * syDis[nr]){
3044dfe5 1015 blst[ir][ic] = kFALSE; continue;
b1957d3c 1016 }
560e5c05 1017 nrow[nr]++; rowId[nr]=ir; kFOUND = kTRUE;
1018 }
1019 if(kFOUND){
1020 vdy[nr].Use(nrow[nr], yres[ir]);
1021 nr++;
b1957d3c 1022 }
b1957d3c 1023 lr = ir; if(nr>=3) break;
29b87567 1024 }
fc0882f3 1025 if(recoParam->GetStreamLevel(AliTRDrecoParam::kTracker) > 3 && fkReconstructor->IsDebugStreaming()){
560e5c05 1026 TTreeSRedirector &cstreamer = *fkReconstructor->GetDebugStream(AliTRDrecoParam::kTracker);
1027 UChar_t stat(0);
1028 if(IsKink()) SETBIT(stat, 1);
1029 if(IsStandAlone()) SETBIT(stat, 2);
1030 cstreamer << "AttachClusters"
1031 << "stat=" << stat
1032 << "det=" << fDet
1033 << "pt=" << fPt
1034 << "s2y=" << s2yTrk
1035 << "r0=" << rowId[0]
1036 << "dy0=" << &vdy[0]
1037 << "m0=" << mean[0]
1038 << "s0=" << syDis[0]
1039 << "r1=" << rowId[1]
1040 << "dy1=" << &vdy[1]
1041 << "m1=" << mean[1]
1042 << "s1=" << syDis[1]
1043 << "r2=" << rowId[2]
1044 << "dy2=" << &vdy[2]
1045 << "m2=" << mean[2]
1046 << "s2=" << syDis[2]
1047 << "\n";
1048 }
1049
1050
1051 // analyze gap in rows attached
1052 if(kRowSelection){
1053 SetErrorMsg(kAttachRowGap);
1054 Int_t rowRemove(-1);
1055 if(nr==2){ // select based on minimum distance to track projection
1056 if(TMath::Abs(mean[0])<TMath::Abs(mean[1])){
1057 if(nrow[1]>nrow[0]) AliDebug(2, Form("Conflicting mean[%f < %f] but ncl[%d < %d].", TMath::Abs(mean[0]), TMath::Abs(mean[1]), nrow[0], nrow[1]));
1058 }else{
1059 if(nrow[1]<nrow[0]) AliDebug(2, Form("Conflicting mean[%f > %f] but ncl[%d > %d].", TMath::Abs(mean[0]), TMath::Abs(mean[1]), nrow[0], nrow[1]));
1060 Swap(nrow[0],nrow[1]); Swap(rowId[0],rowId[1]);
1061 Swap(mean[0],mean[1]); Swap(syDis[0],syDis[1]);
1062 }
1063 rowRemove=1; nr=1;
1064 } else if(nr==3){ // select based on 2 consecutive rows
1065 if(rowId[1]==rowId[0]+1 && rowId[1]!=rowId[2]-1){
1066 nr=2;rowRemove=2;
1067 } else if(rowId[1]!=rowId[0]+1 && rowId[1]==rowId[2]-1){
1068 Swap(nrow[0],nrow[2]); Swap(rowId[0],rowId[2]);
1069 Swap(mean[0],mean[2]); Swap(syDis[0],syDis[2]);
1070 nr=2; rowRemove=2;
1071 }
29b87567 1072 }
560e5c05 1073 if(rowRemove>0){nrow[rowRemove]=0; rowId[rowRemove]=-1;}
29b87567 1074 }
560e5c05 1075 AliDebug(4, Form(" Ncl[%d[%d] + %d[%d] + %d[%d]]", nrow[0], rowId[0], nrow[1], rowId[1], nrow[2], rowId[2]));
1076
1077 if(nr==3){
1078 SetBit(kRowCross, kTRUE); // mark pad row crossing
7c3eecb8 1079 SetErrorMsg(kAttachRow);
560e5c05 1080 const Float_t am[]={TMath::Abs(mean[0]), TMath::Abs(mean[1]), TMath::Abs(mean[2])};
1081 AliDebug(4, Form("complex row configuration\n"
1082 " r[%d] n[%d] m[%6.3f] s[%6.3f]\n"
1083 " r[%d] n[%d] m[%6.3f] s[%6.3f]\n"
1084 " r[%d] n[%d] m[%6.3f] s[%6.3f]\n"
1085 , rowId[0], nrow[0], am[0], syDis[0]
1086 , rowId[1], nrow[1], am[1], syDis[1]
1087 , rowId[2], nrow[2], am[2], syDis[2]));
1088 Int_t id[]={0,1,2}; TMath::Sort(3, am, id, kFALSE);
1089 // backup
1090 Int_t rnn[3]; memcpy(rnn, nrow, 3*sizeof(Int_t));
1091 Int_t rid[3]; memcpy(rid, rowId, 3*sizeof(Int_t));
1092 Double_t rm[3]; memcpy(rm, mean, 3*sizeof(Double_t));
1093 Double_t rs[3]; memcpy(rs, syDis, 3*sizeof(Double_t));
1094 nrow[0]=rnn[id[0]]; rowId[0]=rid[id[0]]; mean[0]=rm[id[0]]; syDis[0]=rs[id[0]];
1095 nrow[1]=rnn[id[1]]; rowId[1]=rid[id[1]]; mean[1]=rm[id[1]]; syDis[1]=rs[id[1]];
1096 nrow[2]=0; rowId[2]=-1; mean[2] = 1.e3; syDis[2] = 1.e3;
1097 AliDebug(4, Form("solved configuration\n"
1098 " r[%d] n[%d] m[%+6.3f] s[%6.3f]\n"
1099 " r[%d] n[%d] m[%+6.3f] s[%6.3f]\n"
1100 " r[%d] n[%d] m[%+6.3f] s[%6.3f]\n"
1101 , rowId[0], nrow[0], mean[0], syDis[0]
1102 , rowId[1], nrow[1], mean[1], syDis[1]
1103 , rowId[2], nrow[2], mean[2], syDis[2]));
1104 nr=2;
1105 } else if(nr==2) {
1106 SetBit(kRowCross, kTRUE); // mark pad row crossing
1107 if(nrow[1] > nrow[0]){ // swap row order
1108 Swap(nrow[0],nrow[1]); Swap(rowId[0],rowId[1]);
1109 Swap(mean[0],mean[1]); Swap(syDis[0],syDis[1]);
1110 }
ee8fb199 1111 }
560e5c05 1112
b1957d3c 1113 // Select and store clusters
1114 // We should consider here :
1115 // 1. How far is the chamber boundary
1116 // 2. How big is the mean
560e5c05 1117 Int_t n(0); Float_t dyc[kNclusters]; memset(dyc,0,kNclusters*sizeof(Float_t));
b1957d3c 1118 for (Int_t ir = 0; ir < nr; ir++) {
560e5c05 1119 Int_t jr(rowId[ir]);
1120 AliDebug(4, Form(" Attaching Ncl[%d]=%d ...", jr, ncl[jr]));
b1957d3c 1121 for (Int_t ic = 0; ic < ncl[jr]; ic++) {
3044dfe5 1122 if(!blst[jr][ic])continue;
1123 c = clst[jr][ic];
560e5c05 1124 Int_t it(c->GetPadTime());
1125 Int_t idx(it+kNtb*ir);
6ad5e6b2 1126 if(fClusters[idx]){
560e5c05 1127 AliDebug(4, Form("Many cluster candidates on row[%2d] tb[%2d].", jr, it));
1128 // TODO should save also the information on where the multiplicity happened and its size
6ad5e6b2 1129 SetErrorMsg(kAttachMultipleCl);
560e5c05 1130 // TODO should also compare with mean and sigma for this row
1131 if(yres[jr][ic] > dyc[idx]) continue;
6ad5e6b2 1132 }
1133
b1957d3c 1134 // TODO proper indexing of clusters !!
6ad5e6b2 1135 fIndexes[idx] = chamber->GetTB(it)->GetGlobalIndex(idxs[jr][ic]);
1136 fClusters[idx] = c;
560e5c05 1137 dyc[idx] = yres[jr][ic];
3e778975 1138 n++;
b1957d3c 1139 }
560e5c05 1140 }
6ad5e6b2 1141 SetN(n);
b1957d3c 1142
29b87567 1143 // number of minimum numbers of clusters expected for the tracklet
6ad5e6b2 1144 if (GetN() < kClmin){
560e5c05 1145 AliDebug(1, Form("NOT ENOUGH CLUSTERS %d ATTACHED TO THE TRACKLET [min %d] FROM FOUND %d.", GetN(), kClmin, n));
7c3eecb8 1146 SetErrorMsg(kAttachClAttach);
e4f2f73d 1147 return kFALSE;
1148 }
0906e73e 1149
e3cf3d02 1150 // Load calibration parameters for this tracklet
1151 Calibrate();
b1957d3c 1152
1153 // calculate dx for time bins in the drift region (calibration aware)
a2abcbc5 1154 Float_t x[2] = {0.,0.}; Int_t tb[2]={0,0};
1155 for (Int_t it = t0, irp=0; irp<2 && it < AliTRDtrackerV1::GetNTimeBins(); it++) {
b1957d3c 1156 if(!fClusters[it]) continue;
1157 x[irp] = fClusters[it]->GetX();
a2abcbc5 1158 tb[irp] = fClusters[it]->GetLocalTimeBin();
b1957d3c 1159 irp++;
e3cf3d02 1160 }
d86ed84c 1161 Int_t dtb = tb[1] - tb[0];
1162 fdX = dtb ? (x[0] - x[1]) / dtb : 0.15;
29b87567 1163 return kTRUE;
e4f2f73d 1164}
1165
03cef9b2 1166//____________________________________________________________
1167void AliTRDseedV1::Bootstrap(const AliTRDReconstructor *rec)
1168{
1169// Fill in all derived information. It has to be called after recovery from file or HLT.
1170// The primitive data are
1171// - list of clusters
1172// - detector (as the detector will be removed from clusters)
1173// - position of anode wire (fX0) - temporary
1174// - track reference position and direction
1175// - momentum of the track
1176// - time bin length [cm]
1177//
1178// A.Bercuci <A.Bercuci@gsi.de> Oct 30th 2008
1179//
4d6aee34 1180 fkReconstructor = rec;
03cef9b2 1181 AliTRDgeometry g;
1182 AliTRDpadPlane *pp = g.GetPadPlane(fDet);
dd8059a8 1183 fPad[0] = pp->GetLengthIPad();
1184 fPad[1] = pp->GetWidthIPad();
525f399b 1185 fPad[2] = TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle());
e3cf3d02 1186 //fSnp = fYref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);
1187 //fTgl = fZref[1];
3e778975 1188 Int_t n = 0, nshare = 0, nused = 0;
03cef9b2 1189 AliTRDcluster **cit = &fClusters[0];
8d2bec9e 1190 for(Int_t ic = kNclusters; ic--; cit++){
03cef9b2 1191 if(!(*cit)) return;
3e778975 1192 n++;
1193 if((*cit)->IsShared()) nshare++;
1194 if((*cit)->IsUsed()) nused++;
03cef9b2 1195 }
3e778975 1196 SetN(n); SetNUsed(nused); SetNShared(nshare);
e3cf3d02 1197 Fit();
03cef9b2 1198 CookLabels();
1199 GetProbability();
1200}
1201
1202
e4f2f73d 1203//____________________________________________________________________
b72f4eaf 1204Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr)
e4f2f73d 1205{
16cca13f 1206//
1207// Linear fit of the clusters attached to the tracklet
1208//
1209// Parameters :
1210// - tilt : switch for tilt pad correction of cluster y position based on
1211// the z, dzdx info from outside [default false].
1212// - zcorr : switch for using z information to correct for anisochronity
1fd9389f 1213// and a finner error parameterization estimation [default false]
16cca13f 1214// Output :
1215// True if successful
1216//
1217// Detailed description
1218//
1219// Fit in the xy plane
1220//
1fd9389f 1221// The fit is performed to estimate the y position of the tracklet and the track
1222// angle in the bending plane. The clusters are represented in the chamber coordinate
1223// system (with respect to the anode wire - see AliTRDtrackerV1::FollowBackProlongation()
1224// on how this is set). The x and y position of the cluster and also their variances
1225// are known from clusterizer level (see AliTRDcluster::GetXloc(), AliTRDcluster::GetYloc(),
1226// AliTRDcluster::GetSX() and AliTRDcluster::GetSY()).
1227// If gaussian approximation is used to calculate y coordinate of the cluster the position
1228// is recalculated taking into account the track angle. The general formula to calculate the
1229// error of cluster position in the gaussian approximation taking into account diffusion and track
1230// inclination is given for TRD by:
1231// BEGIN_LATEX
1232// #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}
1233// END_LATEX
1234//
1235// Since errors are calculated only in the y directions, radial errors (x direction) are mapped to y
1236// by projection i.e.
1237// BEGIN_LATEX
1238// #sigma_{x|y} = tg(#phi) #sigma_{x}
1239// END_LATEX
1240// and also by the lorentz angle correction
1241//
1242// Fit in the xz plane
1243//
1244// The "fit" is performed to estimate the radial position (x direction) where pad row cross happens.
1245// If no pad row crossing the z position is taken from geometry and radial position is taken from the xy
1246// fit (see below).
1247//
1248// There are two methods to estimate the radial position of the pad row cross:
1249// 1. leading cluster radial position : Here the lower part of the tracklet is considered and the last
1250// cluster registered (at radial x0) on this segment is chosen to mark the pad row crossing. The error
1251// of the z estimate is given by :
1252// BEGIN_LATEX
1253// #sigma_{z} = tg(#theta) #Delta x_{x_{0}}/12
1254// END_LATEX
1255// The systematic errors for this estimation are generated by the following sources:
1256// - no charge sharing between pad rows is considered (sharp cross)
1257// - missing cluster at row cross (noise peak-up, under-threshold signal etc.).
1258//
1259// 2. charge fit over the crossing point : Here the full energy deposit along the tracklet is considered
1260// to estimate the position of the crossing by a fit in the qx plane. The errors in the q directions are
1261// parameterized as s_q = q^2. The systematic errors for this estimation are generated by the following sources:
1262// - no general model for the qx dependence
1263// - physical fluctuations of the charge deposit
1264// - gain calibration dependence
1265//
1266// Estimation of the radial position of the tracklet
16cca13f 1267//
1fd9389f 1268// For pad row cross the radial position is taken from the xz fit (see above). Otherwise it is taken as the
1269// interpolation point of the tracklet i.e. the point where the error in y of the fit is minimum. The error
1270// in the y direction of the tracklet is (see AliTRDseedV1::GetCovAt()):
1271// BEGIN_LATEX
1272// #sigma_{y} = #sigma^{2}_{y_{0}} + 2xcov(y_{0}, dy/dx) + #sigma^{2}_{dy/dx}
1273// END_LATEX
1274// and thus the radial position is:
1275// BEGIN_LATEX
1276// x = - cov(y_{0}, dy/dx)/#sigma^{2}_{dy/dx}
1277// END_LATEX
1278//
1279// Estimation of tracklet position error
1280//
1281// The error in y direction is the error of the linear fit at the radial position of the tracklet while in the z
1282// direction is given by the cluster error or pad row cross error. In case of no pad row cross this is given by:
1283// BEGIN_LATEX
1284// #sigma_{y} = #sigma^{2}_{y_{0}} - 2cov^{2}(y_{0}, dy/dx)/#sigma^{2}_{dy/dx} + #sigma^{2}_{dy/dx}
1285// #sigma_{z} = Pad_{length}/12
1286// END_LATEX
1287// For pad row cross the full error is calculated at the radial position of the crossing (see above) and the error
1288// in z by the width of the crossing region - being a matter of parameterization.
1289// BEGIN_LATEX
1290// #sigma_{z} = tg(#theta) #Delta x_{x_{0}}/12
1291// END_LATEX
1292// In case of no tilt correction (default in the barrel tracking) the tilt is taken into account by the rotation of
1293// the covariance matrix. See AliTRDseedV1::GetCovAt() for details.
1294//
1295// Author
1296// A.Bercuci <A.Bercuci@gsi.de>
e4f2f73d 1297
a723055f 1298 if(!fkReconstructor){
1299 AliError("The tracklet needs the reconstruction setup. Please initialize by SetReconstructor().");
1300 return kFALSE;
1301 }
b72f4eaf 1302 if(!IsCalibrated()) Calibrate();
e3cf3d02 1303
29b87567 1304 const Int_t kClmin = 8;
010d62b0 1305
2f7d6ac8 1306 // get track direction
1307 Double_t y0 = fYref[0];
1308 Double_t dydx = fYref[1];
1309 Double_t z0 = fZref[0];
1310 Double_t dzdx = fZref[1];
1311 Double_t yt, zt;
ae4e8b84 1312
5f1ae1e7 1313 AliTRDtrackerV1::AliTRDLeastSquare fitterY;
1314 AliTRDtrackerV1::AliTRDLeastSquare fitterZ;
f301a656 1315
29b87567 1316 // book cluster information
8d2bec9e 1317 Double_t qc[kNclusters], xc[kNclusters], yc[kNclusters], zc[kNclusters], sy[kNclusters];
e3cf3d02 1318
dd8059a8 1319 Int_t n = 0;
4d6aee34 1320 AliTRDcluster *c=NULL, **jc = &fClusters[0];
fc0882f3 1321 const AliTRDrecoParam* const recoParam = fkReconstructor->GetRecoParam(); //the dynamic cast in GetRecoParam is slow, so caching the pointer to it
9eb2d46c 1322 for (Int_t ic=0; ic<kNtb; ic++, ++jc) {
29b87567 1323 xc[ic] = -1.;
1324 yc[ic] = 999.;
1325 zc[ic] = 999.;
1326 sy[ic] = 0.;
9eb2d46c 1327 if(!(c = (*jc))) continue;
29b87567 1328 if(!c->IsInChamber()) continue;
9462866a 1329
29b87567 1330 Float_t w = 1.;
1331 if(c->GetNPads()>4) w = .5;
1332 if(c->GetNPads()>5) w = .2;
010d62b0 1333
1fd9389f 1334 // cluster charge
dd8059a8 1335 qc[n] = TMath::Abs(c->GetQ());
1fd9389f 1336 // pad row of leading
1337
b72f4eaf 1338 // Radial cluster position
e3cf3d02 1339 //Int_t jc = TMath::Max(fN-3, 0);
1340 //xc[fN] = c->GetXloc(fT0, fVD, &qc[jc], &xc[jc]/*, z0 - c->GetX()*dzdx*/);
b72f4eaf 1341 xc[n] = fX0 - c->GetX();
1342
1fd9389f 1343 // extrapolated track to cluster position
dd8059a8 1344 yt = y0 - xc[n]*dydx;
dd8059a8 1345 zt = z0 - xc[n]*dzdx;
1fd9389f 1346
1347 // Recalculate cluster error based on tracking information
1348 c->SetSigmaY2(fS2PRF, fDiffT, fExB, xc[n], zcorr?zt:-1., dydx);
c79857d5 1349 c->SetSigmaZ2(fPad[0]*fPad[0]/12.); // for HLT
1fd9389f 1350 sy[n] = TMath::Sqrt(c->GetSigmaY2());
1351
fc0882f3 1352 yc[n] = recoParam->UseGAUS() ?
1fd9389f 1353 c->GetYloc(y0, sy[n], GetPadWidth()): c->GetY();
1354 zc[n] = c->GetZ();
1355 //optional tilt correction
1356 if(tilt) yc[n] -= (GetTilt()*(zc[n] - zt));
1357
fbe11be7 1358 AliDebug(5, Form(" tb[%2d] dx[%6.3f] y[%6.2f+-%6.3f]", c->GetLocalTimeBin(), xc[n], yc[n], sy[n]));
903326c1 1359 fitterY.AddPoint(&xc[n], yc[n], sy[n]);
0217fcd0 1360 if(IsRowCross()) fitterZ.AddPoint(&xc[n], qc[n], 1.);
dd8059a8 1361 n++;
29b87567 1362 }
3044dfe5 1363
47d5d320 1364 // to few clusters
c79857d5 1365 if (n < kClmin){
1366 SetErrorMsg(kFitFailed);
1367 return kFALSE;
1368 }
2f7d6ac8 1369
d937ad7a 1370 // fit XY
903326c1 1371 if(!fitterY.Eval()){
1372 SetErrorMsg(kFitFailed);
1373 return kFALSE;
1374 }
5f1ae1e7 1375 fYfit[0] = fitterY.GetFunctionParameter(0);
1376 fYfit[1] = -fitterY.GetFunctionParameter(1);
d937ad7a 1377 // store covariance
5f1ae1e7 1378 Double_t p[3];
1379 fitterY.GetCovarianceMatrix(p);
903326c1 1380 fCov[0] = p[1]; // variance of y0
5f1ae1e7 1381 fCov[1] = p[2]; // covariance of y0, dydx
903326c1 1382 fCov[2] = p[0]; // variance of dydx
b1957d3c 1383 // the ref radial position is set at the minimum of
1384 // the y variance of the tracklet
b72f4eaf 1385 fX = -fCov[1]/fCov[2];
903326c1 1386 Float_t xs=fX+.5*AliTRDgeometry::CamHght();
1387 if(xs < 0. || xs > AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght()){
1388 AliDebug(1, Form("Ref radial position ouside chamber x[%5.2f].", fX));
1389 SetErrorMsg(kFitOutside);
1390 return kFALSE;
1391 }
b1957d3c 1392
0217fcd0 1393 // collect second row clusters
1394 Int_t m(0);
b72f4eaf 1395 if(IsRowCross()){
e355f67a 1396/* // THE LEADING CLUSTER METHOD
1fd9389f 1397 Float_t xMin = fX0;
b72f4eaf 1398 Int_t ic=n=kNclusters-1; jc = &fClusters[ic];
1fd9389f 1399 AliTRDcluster *c0 =0x0, **kc = &fClusters[kNtb-1];
1400 for(; ic>kNtb; ic--, --jc, --kc){
1401 if((c0 = (*kc)) && c0->IsInChamber() && (xMin>c0->GetX())) xMin = c0->GetX();
1402 if(!(c = (*jc))) continue;
1403 if(!c->IsInChamber()) continue;
1404 zc[kNclusters-1] = c->GetZ();
1405 fX = fX0 - c->GetX();
1406 }
1407 fZfit[0] = .5*(zc[0]+zc[kNclusters-1]); fZfit[1] = 0.;
1408 // Error parameterization
1409 fS2Z = fdX*fZref[1];
e355f67a 1410 fS2Z *= fS2Z; fS2Z *= 0.2887; // 1/sqrt(12)*/
1411
1fd9389f 1412 // THE FIT X-Q PLANE METHOD
e355f67a 1413 Int_t ic=n=kNclusters-1; jc = &fClusters[ic];
b72f4eaf 1414 for(; ic>kNtb; ic--, --jc){
1415 if(!(c = (*jc))) continue;
1416 if(!c->IsInChamber()) continue;
1417 qc[n] = TMath::Abs(c->GetQ());
1418 xc[n] = fX0 - c->GetX();
1419 zc[n] = c->GetZ();
c79857d5 1420 // Recalculate cluster error based on tracking information
1421 c->SetSigmaY2(fS2PRF, fDiffT, fExB, xc[n], zcorr?(z0 - xc[n]*dzdx):-1., dydx);
1422 c->SetSigmaZ2(fPad[0]*fPad[0]/12.); // for HLT
b72f4eaf 1423 fitterZ.AddPoint(&xc[n], -qc[n], 1.);
0217fcd0 1424 n--;m++;
b72f4eaf 1425 }
0217fcd0 1426 }
1427 // fit XZ
1428 if(m && IsRowCross()){
b72f4eaf 1429 fitterZ.Eval();
5f1ae1e7 1430 if(fitterZ.GetFunctionParameter(1)!=0.){
1431 fX = -fitterZ.GetFunctionParameter(0)/fitterZ.GetFunctionParameter(1);
b72f4eaf 1432 fX=(fX<0.)?0.:fX;
1433 Float_t dl = .5*AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght();
1434 fX=(fX> dl)?dl:fX;
07bbc13c 1435 fX-=.055; // TODO to be understood
b72f4eaf 1436 }
1437
1438 fZfit[0] = .5*(zc[0]+zc[kNclusters-1]); fZfit[1] = 0.;
c850c351 1439 // temporary external error parameterization
1440 fS2Z = 0.05+0.4*TMath::Abs(fZref[1]); fS2Z *= fS2Z;
1441 // TODO correct formula
1442 //fS2Z = sigma_x*TMath::Abs(fZref[1]);
b1957d3c 1443 } else {
0217fcd0 1444 if(IsRowCross() && !m){
1445 AliDebug(1, "Tracklet crossed row but no clusters found in neighbor row.");
1446 }
b1957d3c 1447 fZfit[0] = zc[0]; fZfit[1] = 0.;
dd8059a8 1448 fS2Z = GetPadLength()*GetPadLength()/12.;
29b87567 1449 }
b72f4eaf 1450 fS2Y = fCov[0] +2.*fX*fCov[1] + fX*fX*fCov[2];
29b87567 1451 return kTRUE;
e4f2f73d 1452}
1453
e4f2f73d 1454
f29f13a6 1455/*
e3cf3d02 1456//_____________________________________________________________________________
1457void AliTRDseedV1::FitMI()
1458{
1459//
1460// Fit the seed.
1461// Marian Ivanov's version
1462//
1463// linear fit on the y direction with respect to the reference direction.
1464// The residuals for each x (x = xc - x0) are deduced from:
1465// dy = y - yt (1)
1466// the tilting correction is written :
1467// y = yc + h*(zc-zt) (2)
1468// yt = y0+dy/dx*x (3)
1469// zt = z0+dz/dx*x (4)
1470// from (1),(2),(3) and (4)
1471// dy = yc - y0 - (dy/dx + h*dz/dx)*x + h*(zc-z0)
1472// the last term introduces the correction on y direction due to tilting pads. There are 2 ways to account for this:
1473// 1. use tilting correction for calculating the y
1474// 2. neglect tilting correction here and account for it in the error parametrization of the tracklet.
1475 const Float_t kRatio = 0.8;
1476 const Int_t kClmin = 5;
1477 const Float_t kmaxtan = 2;
1478
1479 if (TMath::Abs(fYref[1]) > kmaxtan){
1480 //printf("Exit: Abs(fYref[1]) = %3.3f, kmaxtan = %3.3f\n", TMath::Abs(fYref[1]), kmaxtan);
1481 return; // Track inclined too much
1482 }
1483
1484 Float_t sigmaexp = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction
dd8059a8 1485 Float_t ycrosscor = GetPadLength() * GetTilt() * 0.5; // Y correction for crossing
e3cf3d02 1486 Int_t fNChange = 0;
1487
1488 Double_t sumw;
1489 Double_t sumwx;
1490 Double_t sumwx2;
1491 Double_t sumwy;
1492 Double_t sumwxy;
1493 Double_t sumwz;
1494 Double_t sumwxz;
1495
1496 // Buffering: Leave it constant fot Performance issues
1497 Int_t zints[kNtb]; // Histograming of the z coordinate
1498 // Get 1 and second max probable coodinates in z
1499 Int_t zouts[2*kNtb];
1500 Float_t allowedz[kNtb]; // Allowed z for given time bin
1501 Float_t yres[kNtb]; // Residuals from reference
dd8059a8 1502 //Float_t anglecor = GetTilt() * fZref[1]; // Correction to the angle
e3cf3d02 1503
1504 Float_t pos[3*kNtb]; memset(pos, 0, 3*kNtb*sizeof(Float_t));
1505 Float_t *fX = &pos[0], *fY = &pos[kNtb], *fZ = &pos[2*kNtb];
1506
1507 Int_t fN = 0; AliTRDcluster *c = 0x0;
1508 fN2 = 0;
1509 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1510 yres[i] = 10000.0;
1511 if (!(c = fClusters[i])) continue;
1512 if(!c->IsInChamber()) continue;
1513 // Residual y
dd8059a8 1514 //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] + GetTilt()*(fZ[i] - fZref[0]);
e3cf3d02 1515 fX[i] = fX0 - c->GetX();
1516 fY[i] = c->GetY();
1517 fZ[i] = c->GetZ();
dd8059a8 1518 yres[i] = fY[i] - GetTilt()*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
e3cf3d02 1519 zints[fN] = Int_t(fZ[i]);
1520 fN++;
1521 }
1522
1523 if (fN < kClmin){
1524 //printf("Exit fN < kClmin: fN = %d\n", fN);
1525 return;
1526 }
1527 Int_t nz = AliTRDtrackerV1::Freq(fN, zints, zouts, kFALSE);
1528 Float_t fZProb = zouts[0];
1529 if (nz <= 1) zouts[3] = 0;
1530 if (zouts[1] + zouts[3] < kClmin) {
1531 //printf("Exit zouts[1] = %d, zouts[3] = %d\n",zouts[1],zouts[3]);
1532 return;
1533 }
1534
1535 // Z distance bigger than pad - length
1536 if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) zouts[3] = 0;
1537
1538 Int_t breaktime = -1;
1539 Bool_t mbefore = kFALSE;
1540 Int_t cumul[kNtb][2];
1541 Int_t counts[2] = { 0, 0 };
1542
1543 if (zouts[3] >= 3) {
1544
1545 //
1546 // Find the break time allowing one chage on pad-rows
1547 // with maximal number of accepted clusters
1548 //
1549 fNChange = 1;
1550 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1551 cumul[i][0] = counts[0];
1552 cumul[i][1] = counts[1];
1553 if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++;
1554 if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++;
1555 }
1556 Int_t maxcount = 0;
1557 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1558 Int_t after = cumul[AliTRDtrackerV1::GetNTimeBins()][0] - cumul[i][0];
1559 Int_t before = cumul[i][1];
1560 if (after + before > maxcount) {
1561 maxcount = after + before;
1562 breaktime = i;
1563 mbefore = kFALSE;
1564 }
1565 after = cumul[AliTRDtrackerV1::GetNTimeBins()-1][1] - cumul[i][1];
1566 before = cumul[i][0];
1567 if (after + before > maxcount) {
1568 maxcount = after + before;
1569 breaktime = i;
1570 mbefore = kTRUE;
1571 }
1572 }
1573 breaktime -= 1;
1574 }
1575
1576 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1577 if (i > breaktime) allowedz[i] = mbefore ? zouts[2] : zouts[0];
1578 if (i <= breaktime) allowedz[i] = (!mbefore) ? zouts[2] : zouts[0];
1579 }
1580
1581 if (((allowedz[0] > allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] < 0)) ||
1582 ((allowedz[0] < allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] > 0))) {
1583 //
1584 // Tracklet z-direction not in correspondance with track z direction
1585 //
1586 fNChange = 0;
1587 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1588 allowedz[i] = zouts[0]; // Only longest taken
1589 }
1590 }
1591
1592 if (fNChange > 0) {
1593 //
1594 // Cross pad -row tracklet - take the step change into account
1595 //
1596 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1597 if (!fClusters[i]) continue;
1598 if(!fClusters[i]->IsInChamber()) continue;
1599 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
1600 // Residual y
dd8059a8 1601 //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] + GetTilt()*(fZ[i] - fZref[0]);
1602 yres[i] = fY[i] - GetTilt()*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
f29f13a6 1603// if (TMath::Abs(fZ[i] - fZProb) > 2) {
dd8059a8 1604// if (fZ[i] > fZProb) yres[i] += GetTilt() * GetPadLength();
1605// if (fZ[i] < fZProb) yres[i] -= GetTilt() * GetPadLength();
f29f13a6 1606 }
e3cf3d02 1607 }
1608 }
1609
1610 Double_t yres2[kNtb];
1611 Double_t mean;
1612 Double_t sigma;
1613 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1614 if (!fClusters[i]) continue;
1615 if(!fClusters[i]->IsInChamber()) continue;
1616 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
1617 yres2[fN2] = yres[i];
1618 fN2++;
1619 }
1620 if (fN2 < kClmin) {
1621 //printf("Exit fN2 < kClmin: fN2 = %d\n", fN2);
1622 fN2 = 0;
1623 return;
1624 }
1625 AliMathBase::EvaluateUni(fN2,yres2,mean,sigma, Int_t(fN2*kRatio-2.));
1626 if (sigma < sigmaexp * 0.8) {
1627 sigma = sigmaexp;
1628 }
1629 //Float_t fSigmaY = sigma;
1630
1631 // Reset sums
1632 sumw = 0;
1633 sumwx = 0;
1634 sumwx2 = 0;
1635 sumwy = 0;
1636 sumwxy = 0;
1637 sumwz = 0;
1638 sumwxz = 0;
1639
1640 fN2 = 0;
1641 Float_t fMeanz = 0;
1642 Float_t fMPads = 0;
1643 fUsable = 0;
1644 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1645 if (!fClusters[i]) continue;
1646 if (!fClusters[i]->IsInChamber()) continue;
1647 if (TMath::Abs(fZ[i] - allowedz[i]) > 2){fClusters[i] = 0x0; continue;}
1648 if (TMath::Abs(yres[i] - mean) > 4.0 * sigma){fClusters[i] = 0x0; continue;}
1649 SETBIT(fUsable,i);
1650 fN2++;
1651 fMPads += fClusters[i]->GetNPads();
1652 Float_t weight = 1.0;
1653 if (fClusters[i]->GetNPads() > 4) weight = 0.5;
1654 if (fClusters[i]->GetNPads() > 5) weight = 0.2;
1655
1656
1657 Double_t x = fX[i];
1658 //printf("x = %7.3f dy = %7.3f fit %7.3f\n", x, yres[i], fY[i]-yres[i]);
1659
1660 sumw += weight;
1661 sumwx += x * weight;
1662 sumwx2 += x*x * weight;
1663 sumwy += weight * yres[i];
1664 sumwxy += weight * (yres[i]) * x;
1665 sumwz += weight * fZ[i];
1666 sumwxz += weight * fZ[i] * x;
1667
1668 }
1669
1670 if (fN2 < kClmin){
1671 //printf("Exit fN2 < kClmin(2): fN2 = %d\n",fN2);
1672 fN2 = 0;
1673 return;
1674 }
1675 fMeanz = sumwz / sumw;
1676 Float_t correction = 0;
1677 if (fNChange > 0) {
1678 // Tracklet on boundary
1679 if (fMeanz < fZProb) correction = ycrosscor;
1680 if (fMeanz > fZProb) correction = -ycrosscor;
1681 }
1682
1683 Double_t det = sumw * sumwx2 - sumwx * sumwx;
1684 fYfit[0] = (sumwx2 * sumwy - sumwx * sumwxy) / det;
1685 fYfit[1] = (sumw * sumwxy - sumwx * sumwy) / det;
1686
1687 fS2Y = 0;
1688 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1689 if (!TESTBIT(fUsable,i)) continue;
1690 Float_t delta = yres[i] - fYfit[0] - fYfit[1] * fX[i];
1691 fS2Y += delta*delta;
1692 }
1693 fS2Y = TMath::Sqrt(fS2Y / Float_t(fN2-2));
1694 // TEMPORARY UNTIL covariance properly calculated
1695 fS2Y = TMath::Max(fS2Y, Float_t(.1));
1696
1697 fZfit[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det;
1698 fZfit[1] = (sumw * sumwxz - sumwx * sumwz) / det;
1699// fYfitR[0] += fYref[0] + correction;
1700// fYfitR[1] += fYref[1];
1701// fYfit[0] = fYfitR[0];
1702 fYfit[1] = -fYfit[1];
1703
1704 UpdateUsed();
f29f13a6 1705}*/
e3cf3d02 1706
e4f2f73d 1707//___________________________________________________________________
203967fc 1708void AliTRDseedV1::Print(Option_t *o) const
e4f2f73d 1709{
1710 //
1711 // Printing the seedstatus
1712 //
1713
b72f4eaf 1714 AliInfo(Form("Det[%3d] X0[%7.2f] Pad{L[%5.2f] W[%5.2f] Tilt[%+6.2f]}", fDet, fX0, GetPadLength(), GetPadWidth(), GetTilt()));
dd8059a8 1715 AliInfo(Form("N[%2d] Nused[%2d] Nshared[%2d] [%d]", GetN(), GetNUsed(), GetNShared(), fN));
b72f4eaf 1716 AliInfo(Form("FLAGS : RC[%c] Kink[%c] SA[%c]", IsRowCross()?'y':'n', IsKink()?'y':'n', IsStandAlone()?'y':'n'));
525f399b 1717 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 1718
1719 Double_t cov[3], x=GetX();
1720 GetCovAt(x, cov);
1721 AliInfo(" | x[cm] | y[cm] | z[cm] | dydx | dzdx |");
1722 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]));
16cca13f 1723 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 1724 AliInfo(Form("P / Pt [GeV/c] = %f / %f", GetMomentum(), fPt));
68f9b6bd 1725 if(IsStandAlone()) AliInfo(Form("C Rieman / Vertex [1/cm] = %f / %f", fC[0], fC[1]));
ee8fb199 1726 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]));
1727 AliInfo(Form("PID = %5.3f / %5.3f / %5.3f / %5.3f / %5.3f", fProb[0], fProb[1], fProb[2], fProb[3], fProb[4]));
203967fc 1728
1729 if(strcmp(o, "a")!=0) return;
1730
4dc4dc2e 1731 AliTRDcluster* const* jc = &fClusters[0];
8d2bec9e 1732 for(int ic=0; ic<kNclusters; ic++, jc++) {
4dc4dc2e 1733 if(!(*jc)) continue;
203967fc 1734 (*jc)->Print(o);
4dc4dc2e 1735 }
e4f2f73d 1736}
47d5d320 1737
203967fc 1738
1739//___________________________________________________________________
1740Bool_t AliTRDseedV1::IsEqual(const TObject *o) const
1741{
1742 // Checks if current instance of the class has the same essential members
1743 // as the given one
1744
1745 if(!o) return kFALSE;
1746 const AliTRDseedV1 *inTracklet = dynamic_cast<const AliTRDseedV1*>(o);
1747 if(!inTracklet) return kFALSE;
1748
1749 for (Int_t i = 0; i < 2; i++){
e3cf3d02 1750 if ( fYref[i] != inTracklet->fYref[i] ) return kFALSE;
1751 if ( fZref[i] != inTracklet->fZref[i] ) return kFALSE;
203967fc 1752 }
1753
e3cf3d02 1754 if ( fS2Y != inTracklet->fS2Y ) return kFALSE;
dd8059a8 1755 if ( GetTilt() != inTracklet->GetTilt() ) return kFALSE;
1756 if ( GetPadLength() != inTracklet->GetPadLength() ) return kFALSE;
203967fc 1757
8d2bec9e 1758 for (Int_t i = 0; i < kNclusters; i++){
e3cf3d02 1759// if ( fX[i] != inTracklet->GetX(i) ) return kFALSE;
1760// if ( fY[i] != inTracklet->GetY(i) ) return kFALSE;
1761// if ( fZ[i] != inTracklet->GetZ(i) ) return kFALSE;
1762 if ( fIndexes[i] != inTracklet->fIndexes[i] ) return kFALSE;
203967fc 1763 }
f29f13a6 1764// if ( fUsable != inTracklet->fUsable ) return kFALSE;
203967fc 1765
1766 for (Int_t i=0; i < 2; i++){
e3cf3d02 1767 if ( fYfit[i] != inTracklet->fYfit[i] ) return kFALSE;
1768 if ( fZfit[i] != inTracklet->fZfit[i] ) return kFALSE;
1769 if ( fLabels[i] != inTracklet->fLabels[i] ) return kFALSE;
203967fc 1770 }
1771
e3cf3d02 1772/* if ( fMeanz != inTracklet->GetMeanz() ) return kFALSE;
1773 if ( fZProb != inTracklet->GetZProb() ) return kFALSE;*/
3e778975 1774 if ( fN != inTracklet->fN ) return kFALSE;
1775 //if ( fNUsed != inTracklet->fNUsed ) return kFALSE;
e3cf3d02 1776 //if ( fFreq != inTracklet->GetFreq() ) return kFALSE;
1777 //if ( fNChange != inTracklet->GetNChange() ) return kFALSE;
203967fc 1778
e3cf3d02 1779 if ( fC != inTracklet->fC ) return kFALSE;
1780 //if ( fCC != inTracklet->GetCC() ) return kFALSE;
1781 if ( fChi2 != inTracklet->fChi2 ) return kFALSE;
203967fc 1782 // if ( fChi2Z != inTracklet->GetChi2Z() ) return kFALSE;
1783
e3cf3d02 1784 if ( fDet != inTracklet->fDet ) return kFALSE;
b25a5e9e 1785 if ( fPt != inTracklet->fPt ) return kFALSE;
e3cf3d02 1786 if ( fdX != inTracklet->fdX ) return kFALSE;
203967fc 1787
8d2bec9e 1788 for (Int_t iCluster = 0; iCluster < kNclusters; iCluster++){
203967fc 1789 AliTRDcluster *curCluster = fClusters[iCluster];
e3cf3d02 1790 AliTRDcluster *inCluster = inTracklet->fClusters[iCluster];
203967fc 1791 if (curCluster && inCluster){
1792 if (! curCluster->IsEqual(inCluster) ) {
1793 curCluster->Print();
1794 inCluster->Print();
1795 return kFALSE;
1796 }
1797 } else {
1798 // if one cluster exists, and corresponding
1799 // in other tracklet doesn't - return kFALSE
1800 if(curCluster || inCluster) return kFALSE;
1801 }
1802 }
1803 return kTRUE;
1804}
5d401b45 1805