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