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