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