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