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