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