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