changes AliTRDSignalIndex
[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)
0906e73e 80 ,fMom(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
6e4d4425 105 fRefCov[0] = 1.; fRefCov[1] = 0.; fRefCov[2] = 1.;
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)
e3cf3d02 126 ,fMom(0.)
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;
29b87567 203 target.fMom = fMom;
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));
224 memcpy(target.fRefCov, fRefCov, 3*sizeof(Double_t));
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;
b1957d3c 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;
e3cf3d02 266 fMom=0.;
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
282 fRefCov[0] = 1.; fRefCov[1] = 0.; fRefCov[2] = 1.;
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//____________________________________________________________________
b1957d3c 289void AliTRDseedV1::UpDate(const AliTRDtrackV1 *trk)
290{
291 // update tracklet reference position from the TRD track
292 // Funny name to avoid the clash with the function AliTRDseed::Update() (has to be made obselete)
293
e3cf3d02 294 Double_t fSnp = trk->GetSnp();
295 Double_t fTgl = trk->GetTgl();
b1957d3c 296 fMom = trk->GetP();
297 fYref[1] = fSnp/(1. - fSnp*fSnp);
298 fZref[1] = fTgl;
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//____________________________________________________________________
454void AliTRDseedV1::GetClusterXY(const AliTRDcluster *c, Double_t &x, Double_t &y)
455{
456// Return corrected position of the cluster taking into
457// account variation of the drift velocity with drift length.
458
459
460 // drift velocity correction TODO to be moved to the clusterizer
461 const Float_t cx[] = {
462 -9.6280e-02, 1.3091e-01,-1.7415e-02,-9.9221e-02,-1.2040e-01,-9.5493e-02,
463 -5.0041e-02,-1.6726e-02, 3.5756e-03, 1.8611e-02, 2.6378e-02, 3.3823e-02,
464 3.4811e-02, 3.5282e-02, 3.5386e-02, 3.6047e-02, 3.5201e-02, 3.4384e-02,
465 3.2864e-02, 3.1932e-02, 3.2051e-02, 2.2539e-02,-2.5154e-02,-1.2050e-01,
466 -1.2050e-01
467 };
468
469 // PRF correction TODO to be replaced by the gaussian
470 // approximation with full error parametrization and // moved to the clusterizer
471 const Float_t cy[AliTRDgeometry::kNlayer][3] = {
472 { 4.014e-04, 8.605e-03, -6.880e+00},
473 {-3.061e-04, 9.663e-03, -6.789e+00},
474 { 1.124e-03, 1.105e-02, -6.825e+00},
475 {-1.527e-03, 1.231e-02, -6.777e+00},
476 { 2.150e-03, 1.387e-02, -6.783e+00},
477 {-1.296e-03, 1.486e-02, -6.825e+00}
478 };
479
480 Int_t ily = AliTRDgeometry::GetLayer(c->GetDetector());
481 x = c->GetX() - cx[c->GetLocalTimeBin()];
482 y = c->GetY() + cy[ily][0] + cy[ily][1] * TMath::Sin(cy[ily][2] * c->GetCenter());
483 return;
484}
b83573da 485
bcb6fb78 486//____________________________________________________________________
487Float_t AliTRDseedV1::GetdQdl(Int_t ic) const
488{
3ee48d6e 489// Using the linear approximation of the track inside one TRD chamber (TRD tracklet)
490// the charge per unit length can be written as:
491// BEGIN_LATEX
492// #frac{dq}{dl} = #frac{q_{c}}{dx * #sqrt{1 + #(){#frac{dy}{dx}}^{2}_{fit} + #(){#frac{dy}{dx}}^{2}_{ref}}}
493// END_LATEX
494// where qc is the total charge collected in the current time bin and dx is the length
495// of the time bin. For the moment (Jan 20 2009) only pad row cross corrections are
496// considered for the charge but none are applied for drift velocity variations along
497// the drift region or assymetry of the TRF
498//
499// Author : Alex Bercuci <A.Bercuci@gsi.de>
500//
501 Float_t dq = 0.;
502 if(fClusters[ic]) dq += TMath::Abs(fClusters[ic]->GetQ());
8e709c82 503 if(fClusters[ic+kNtb]) dq += TMath::Abs(fClusters[ic+kNtb]->GetQ());
504 if(dq<1.e-3 || fdX < 1.e-3) return 0.;
3ee48d6e 505
506 return dq/fdX/TMath::Sqrt(1. + fYfit[1]*fYfit[1] + fZref[1]*fZref[1]);
bcb6fb78 507}
508
0906e73e 509//____________________________________________________________________
3e778975 510Float_t* AliTRDseedV1::GetProbability(Bool_t force)
0906e73e 511{
3e778975 512 if(!force) return &fProb[0];
513 if(!CookPID()) return 0x0;
514 return &fProb[0];
515}
516
517//____________________________________________________________
518Bool_t AliTRDseedV1::CookPID()
519{
0906e73e 520// Fill probability array for tracklet from the DB.
521//
522// Parameters
523//
524// Output
525// returns pointer to the probability array and 0x0 if missing DB access
526//
527// Detailed description
528
29b87567 529
530 // retrive calibration db
0906e73e 531 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
532 if (!calibration) {
533 AliError("No access to calibration data");
3e778975 534 return kFALSE;
0906e73e 535 }
536
3a039a31 537 if (!fReconstructor) {
538 AliError("Reconstructor not set.");
3e778975 539 return kFALSE;
4ba1d6ae 540 }
541
0906e73e 542 // Retrieve the CDB container class with the parametric detector response
3a039a31 543 const AliTRDCalPID *pd = calibration->GetPIDObject(fReconstructor->GetPIDMethod());
0906e73e 544 if (!pd) {
545 AliError("No access to AliTRDCalPID object");
3e778975 546 return kFALSE;
0906e73e 547 }
29b87567 548 //AliInfo(Form("Method[%d] : %s", fReconstructor->GetRecoParam() ->GetPIDMethod(), pd->IsA()->GetName()));
10f75631 549
29b87567 550 // calculate tracklet length TO DO
0906e73e 551 Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
552 /// TMath::Sqrt((1.0 - fSnp[iPlane]*fSnp[iPlane]) / (1.0 + fTgl[iPlane]*fTgl[iPlane]));
553
554 //calculate dE/dx
3a039a31 555 CookdEdx(fReconstructor->GetNdEdxSlices());
0906e73e 556
557 // Sets the a priori probabilities
558 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
ae4e8b84 559 fProb[ispec] = pd->GetProbability(ispec, fMom, &fdEdx[0], length, GetPlane());
0906e73e 560 }
561
3e778975 562 return kTRUE;
0906e73e 563}
564
565//____________________________________________________________________
e4f2f73d 566Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
567{
568 //
569 // Returns a quality measurement of the current seed
570 //
571
dd8059a8 572 Float_t zcorr = kZcorr ? GetTilt() * (fZfit[0] - fZref[0]) : 0.;
29b87567 573 return
3e778975 574 .5 * TMath::Abs(18.0 - GetN())
29b87567 575 + 10.* TMath::Abs(fYfit[1] - fYref[1])
576 + 5. * TMath::Abs(fYfit[0] - fYref[0] + zcorr)
dd8059a8 577 + 2. * TMath::Abs(fZfit[0] - fZref[0]) / GetPadLength();
e4f2f73d 578}
579
580//____________________________________________________________________
d937ad7a 581void AliTRDseedV1::GetCovAt(Double_t x, Double_t *cov) const
0906e73e 582{
d937ad7a 583// Computes covariance in the y-z plane at radial point x (in tracking coordinates)
584// and returns the results in the preallocated array cov[3] as :
585// cov[0] = Var(y)
586// cov[1] = Cov(yz)
587// cov[2] = Var(z)
588//
589// Details
590//
591// For the linear transformation
592// BEGIN_LATEX
593// Y = T_{x} X^{T}
594// END_LATEX
595// The error propagation has the general form
596// BEGIN_LATEX
597// C_{Y} = T_{x} C_{X} T_{x}^{T}
598// END_LATEX
599// We apply this formula 2 times. First to calculate the covariance of the tracklet
600// at point x we consider:
601// BEGIN_LATEX
602// T_{x} = (1 x); X=(y0 dy/dx); C_{X}=#(){#splitline{Var(y0) Cov(y0, dy/dx)}{Cov(y0, dy/dx) Var(dy/dx)}}
603// END_LATEX
604// and secondly to take into account the tilt angle
605// BEGIN_LATEX
606// T_{#alpha} = #(){#splitline{cos(#alpha) __ sin(#alpha)}{-sin(#alpha) __ cos(#alpha)}}; X=(y z); C_{X}=#(){#splitline{Var(y) 0}{0 Var(z)}}
607// END_LATEX
608//
609// using simple trigonometrics one can write for this last case
610// BEGIN_LATEX
611// 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})}}
612// END_LATEX
613// which can be aproximated for small alphas (2 deg) with
614// BEGIN_LATEX
615// 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}}}
616// END_LATEX
617//
618// before applying the tilt rotation we also apply systematic uncertainties to the tracklet
619// position which can be tunned from outside via the AliTRDrecoParam::SetSysCovMatrix(). They might
620// account for extra misalignment/miscalibration uncertainties.
621//
622// Author :
623// Alex Bercuci <A.Bercuci@gsi.de>
624// Date : Jan 8th 2009
625//
b1957d3c 626
627
d937ad7a 628 Double_t xr = fX0-x;
629 Double_t sy2 = fCov[0] +2.*xr*fCov[1] + xr*xr*fCov[2];
dd8059a8 630 Double_t sz2 = GetPadLength()*GetPadLength()/12.;
0906e73e 631
d937ad7a 632 // insert systematic uncertainties
633 Double_t sys[15];
634 fReconstructor->GetRecoParam()->GetSysCovMatrix(sys);
635 sy2 += sys[0];
636 sz2 += sys[1];
637
638 // rotate covariance matrix
dd8059a8 639 Double_t t2 = GetTilt()*GetTilt();
d937ad7a 640 Double_t correction = 1./(1. + t2);
641 cov[0] = (sy2+t2*sz2)*correction;
dd8059a8 642 cov[1] = GetTilt()*(sz2 - sy2)*correction;
d937ad7a 643 cov[2] = (t2*sy2+sz2)*correction;
644}
eb38ed55 645
0906e73e 646
d937ad7a 647//____________________________________________________________________
e3cf3d02 648void AliTRDseedV1::Calibrate()
d937ad7a 649{
e3cf3d02 650// Retrieve calibration and position parameters from OCDB.
651// The following information are used
d937ad7a 652// - detector index
e3cf3d02 653// - column and row position of first attached cluster. If no clusters are attached
654// to the tracklet a random central chamber position (c=70, r=7) will be used.
655//
656// The following information is cached in the tracklet
657// t0 (trigger delay)
658// drift velocity
659// PRF width
660// omega*tau = tg(a_L)
661// diffusion coefficients (longitudinal and transversal)
d937ad7a 662//
663// Author :
664// Alex Bercuci <A.Bercuci@gsi.de>
665// Date : Jan 8th 2009
666//
eb38ed55 667
d937ad7a 668 AliCDBManager *cdb = AliCDBManager::Instance();
669 if(cdb->GetRun() < 0){
670 AliError("OCDB manager not properly initialized");
671 return;
672 }
0906e73e 673
e3cf3d02 674 AliTRDcalibDB *calib = AliTRDcalibDB::Instance();
675 AliTRDCalROC *vdROC = calib->GetVdriftROC(fDet),
676 *t0ROC = calib->GetT0ROC(fDet);;
677 const AliTRDCalDet *vdDet = calib->GetVdriftDet();
678 const AliTRDCalDet *t0Det = calib->GetT0Det();
d937ad7a 679
680 Int_t col = 70, row = 7;
681 AliTRDcluster **c = &fClusters[0];
3e778975 682 if(GetN()){
d937ad7a 683 Int_t ic = 0;
8d2bec9e 684 while (ic<kNclusters && !(*c)){ic++; c++;}
d937ad7a 685 if(*c){
686 col = (*c)->GetPadCol();
687 row = (*c)->GetPadRow();
688 }
689 }
3a039a31 690
e3cf3d02 691 fT0 = t0Det->GetValue(fDet) + t0ROC->GetValue(col,row);
692 fVD = vdDet->GetValue(fDet) * vdROC->GetValue(col, row);
693 fS2PRF = calib->GetPRFWidth(fDet, col, row); fS2PRF *= fS2PRF;
694 fExB = AliTRDCommonParam::Instance()->GetOmegaTau(fVD);
695 AliTRDCommonParam::Instance()->GetDiffCoeff(fDiffL,
696 fDiffT, fVD);
697 SetBit(kCalib, kTRUE);
0906e73e 698}
699
0906e73e 700//____________________________________________________________________
29b87567 701void AliTRDseedV1::SetOwner()
0906e73e 702{
29b87567 703 //AliInfo(Form("own [%s] fOwner[%s]", own?"YES":"NO", fOwner?"YES":"NO"));
704
705 if(TestBit(kOwner)) return;
8d2bec9e 706 for(int ic=0; ic<kNclusters; ic++){
29b87567 707 if(!fClusters[ic]) continue;
708 fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
709 }
710 SetBit(kOwner);
0906e73e 711}
712
eb2b4f91 713//____________________________________________________________
714void AliTRDseedV1::SetPadPlane(AliTRDpadPlane *p)
715{
716// Shortcut method to initialize pad geometry.
717 if(!p) return;
718 SetTilt(TMath::Tan(TMath::DegToRad()*p->GetTiltingAngle()));
719 SetPadLength(p->GetLengthIPad());
720 SetPadWidth(p->GetWidthIPad());
721}
722
723
f29f13a6 724// //____________________________________________________________________
725// Bool_t AliTRDseedV1::AttachClustersIter(AliTRDtrackingChamber *chamber, Float_t quality, Bool_t kZcorr, AliTRDcluster *c)
726// {
727// //
728// // Iterative process to register clusters to the seed.
729// // In iteration 0 we try only one pad-row and if quality not
730// // sufficient we try 2 pad-rows (about 5% of tracks cross 2 pad-rows)
731// //
732// // debug level 7
733// //
734//
735// if(!fReconstructor->GetRecoParam() ){
736// AliError("Seed can not be used without a valid RecoParam.");
737// return kFALSE;
738// }
739//
740// AliTRDchamberTimeBin *layer = 0x0;
741// if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7){
742// AliTRDtrackingChamber ch(*chamber);
743// ch.SetOwner();
744// TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker);
745// cstreamer << "AttachClustersIter"
746// << "chamber.=" << &ch
747// << "tracklet.=" << this
748// << "\n";
749// }
750//
751// Float_t tquality;
752// Double_t kroady = fReconstructor->GetRecoParam() ->GetRoad1y();
dd8059a8 753// Double_t kroadz = GetPadLength() * .5 + 1.;
f29f13a6 754//
755// // initialize configuration parameters
dd8059a8 756// Float_t zcorr = kZcorr ? GetTilt() * (fZfit[0] - fZref[0]) : 0.;
f29f13a6 757// Int_t niter = kZcorr ? 1 : 2;
758//
759// Double_t yexp, zexp;
760// Int_t ncl = 0;
761// // start seed update
762// for (Int_t iter = 0; iter < niter; iter++) {
763// ncl = 0;
764// for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
765// if(!(layer = chamber->GetTB(iTime))) continue;
766// if(!Int_t(*layer)) continue;
767//
768// // define searching configuration
769// Double_t dxlayer = layer->GetX() - fX0;
770// if(c){
771// zexp = c->GetZ();
772// //Try 2 pad-rows in second iteration
773// if (iter > 0) {
774// zexp = fZref[0] + fZref[1] * dxlayer - zcorr;
dd8059a8 775// if (zexp > c->GetZ()) zexp = c->GetZ() + GetPadLength()*0.5;
776// if (zexp < c->GetZ()) zexp = c->GetZ() - GetPadLength()*0.5;
f29f13a6 777// }
778// } else zexp = fZref[0] + (kZcorr ? fZref[1] * dxlayer : 0.);
779// yexp = fYref[0] + fYref[1] * dxlayer - zcorr;
780//
781// // Get and register cluster
782// Int_t index = layer->SearchNearestCluster(yexp, zexp, kroady, kroadz);
783// if (index < 0) continue;
784// AliTRDcluster *cl = (*layer)[index];
785//
786// fIndexes[iTime] = layer->GetGlobalIndex(index);
787// fClusters[iTime] = cl;
788// // fY[iTime] = cl->GetY();
789// // fZ[iTime] = cl->GetZ();
790// ncl++;
791// }
792// if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d ncl [%d] = %d", iter, fDet, ncl));
793//
794// if(ncl>1){
795// // calculate length of the time bin (calibration aware)
796// Int_t irp = 0; Float_t x[2]={0., 0.}; Int_t tb[2] = {0,0};
e3cf3d02 797// for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
798// if(!fClusters[iTime]) continue;
f29f13a6 799// x[irp] = fClusters[iTime]->GetX();
800// tb[irp] = iTime;
801// irp++;
802// if(irp==2) break;
e3cf3d02 803// }
f29f13a6 804// Int_t dtb = tb[1] - tb[0];
805// fdX = dtb ? (x[0] - x[1]) / dtb : 0.15;
806//
807// // update X0 from the clusters (calibration/alignment aware)
808// for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
809// if(!(layer = chamber->GetTB(iTime))) continue;
810// if(!layer->IsT0()) continue;
811// if(fClusters[iTime]){
812// fX0 = fClusters[iTime]->GetX();
813// break;
814// } else { // we have to infere the position of the anode wire from the other clusters
815// for (Int_t jTime = iTime+1; jTime < AliTRDtrackerV1::GetNTimeBins(); jTime++) {
816// if(!fClusters[jTime]) continue;
817// fX0 = fClusters[jTime]->GetX() + fdX * (jTime - iTime);
818// break;
819// }
820// }
821// }
822//
823// // update YZ reference point
824// // TODO
825//
826// // update x reference positions (calibration/alignment aware)
827// // for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
828// // if(!fClusters[iTime]) continue;
829// // fX[iTime] = fX0 - fClusters[iTime]->GetX();
830// // }
831//
832// FitMI();
833// }
834// if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d nclFit [%d] = %d", iter, fDet, fN2));
835//
836// if(IsOK()){
837// tquality = GetQuality(kZcorr);
838// if(tquality < quality) break;
839// else quality = tquality;
840// }
841// kroadz *= 2.;
842// } // Loop: iter
843// if (!IsOK()) return kFALSE;
844//
845// if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=1) CookLabels();
846//
847// // load calibration params
848// Calibrate();
849// UpdateUsed();
850// return kTRUE;
851// }
e4f2f73d 852
853//____________________________________________________________________
b1957d3c 854Bool_t AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *chamber, Bool_t tilt)
e4f2f73d 855{
856 //
857 // Projective algorithm to attach clusters to seeding tracklets
858 //
859 // Parameters
860 //
861 // Output
862 //
863 // Detailed description
864 // 1. Collapse x coordinate for the full detector plane
865 // 2. truncated mean on y (r-phi) direction
866 // 3. purge clusters
867 // 4. truncated mean on z direction
868 // 5. purge clusters
869 // 6. fit tracklet
870 //
b1957d3c 871 Bool_t kPRINT = kFALSE;
29b87567 872 if(!fReconstructor->GetRecoParam() ){
873 AliError("Seed can not be used without a valid RecoParam.");
874 return kFALSE;
875 }
b1957d3c 876 // Initialize reco params for this tracklet
877 // 1. first time bin in the drift region
878 Int_t t0 = 4;
879 Int_t kClmin = Int_t(fReconstructor->GetRecoParam() ->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
29b87567 880
b1957d3c 881 Double_t syRef = TMath::Sqrt(fRefCov[0]);
29b87567 882 //define roads
b1957d3c 883 Double_t kroady = 1.;
884 //fReconstructor->GetRecoParam() ->GetRoad1y();
dd8059a8 885 Double_t kroadz = GetPadLength() * 1.5 + 1.;
b1957d3c 886 if(kPRINT) printf("AttachClusters() sy[%f] road[%f]\n", syRef, kroady);
29b87567 887
888 // working variables
b1957d3c 889 const Int_t kNrows = 16;
8d2bec9e 890 AliTRDcluster *clst[kNrows][kNclusters];
b1957d3c 891 Double_t cond[4], dx, dy, yt, zt,
8d2bec9e 892 yres[kNrows][kNclusters];
893 Int_t idxs[kNrows][kNclusters], ncl[kNrows], ncls = 0;
b1957d3c 894 memset(ncl, 0, kNrows*sizeof(Int_t));
8d2bec9e 895 memset(clst, 0, kNrows*kNclusters*sizeof(AliTRDcluster*));
b1957d3c 896
29b87567 897 // Do cluster projection
b1957d3c 898 AliTRDcluster *c = 0x0;
29b87567 899 AliTRDchamberTimeBin *layer = 0x0;
b1957d3c 900 Bool_t kBUFFER = kFALSE;
901 for (Int_t it = 0; it < AliTRDtrackerV1::GetNTimeBins(); it++) {
902 if(!(layer = chamber->GetTB(it))) continue;
29b87567 903 if(!Int_t(*layer)) continue;
904
b1957d3c 905 dx = fX0 - layer->GetX();
906 yt = fYref[0] - fYref[1] * dx;
907 zt = fZref[0] - fZref[1] * dx;
908 if(kPRINT) printf("\t%2d dx[%f] yt[%f] zt[%f]\n", it, dx, yt, zt);
909
910 // select clusters on a 5 sigmaKalman level
911 cond[0] = yt; cond[2] = kroady;
912 cond[1] = zt; cond[3] = kroadz;
913 Int_t n=0, idx[6];
914 layer->GetClusters(cond, idx, n, 6);
915 for(Int_t ic = n; ic--;){
916 c = (*layer)[idx[ic]];
917 dy = yt - c->GetY();
dd8059a8 918 dy += tilt ? GetTilt() * (c->GetZ() - zt) : 0.;
b1957d3c 919 // select clusters on a 3 sigmaKalman level
920/* if(tilt && TMath::Abs(dy) > 3.*syRef){
921 printf("too large !!!\n");
922 continue;
923 }*/
924 Int_t r = c->GetPadRow();
925 if(kPRINT) printf("\t\t%d dy[%f] yc[%f] r[%d]\n", ic, TMath::Abs(dy), c->GetY(), r);
926 clst[r][ncl[r]] = c;
927 idxs[r][ncl[r]] = idx[ic];
928 yres[r][ncl[r]] = dy;
929 ncl[r]++; ncls++;
930
8d2bec9e 931 if(ncl[r] >= kNclusters) {
932 AliWarning(Form("Cluster candidates reached limit %d. Some may be lost.", kNclusters));
b1957d3c 933 kBUFFER = kTRUE;
29b87567 934 break;
935 }
936 }
b1957d3c 937 if(kBUFFER) break;
29b87567 938 }
b1957d3c 939 if(kPRINT) printf("Found %d clusters\n", ncls);
940 if(ncls<kClmin) return kFALSE;
941
942 // analyze each row individualy
943 Double_t mean, syDis;
944 Int_t nrow[] = {0, 0, 0}, nr = 0, lr=-1;
945 for(Int_t ir=kNrows; ir--;){
946 if(!(ncl[ir])) continue;
947 if(lr>0 && lr-ir != 1){
948 if(kPRINT) printf("W - gap in rows attached !!\n");
29b87567 949 }
b1957d3c 950 if(kPRINT) printf("\tir[%d] lr[%d] n[%d]\n", ir, lr, ncl[ir]);
951 // Evaluate truncated mean on the y direction
952 if(ncl[ir] > 3) AliMathBase::EvaluateUni(ncl[ir], yres[ir], mean, syDis, Int_t(ncl[ir]*.8));
953 else {
954 mean = 0.; syDis = 0.;
955 }
956
957 // TODO check mean and sigma agains cluster resolution !!
958 if(kPRINT) printf("\tr[%2d] m[%f %5.3fsigma] s[%f]\n", ir, mean, TMath::Abs(mean/syRef), syDis);
959 // select clusters on a 3 sigmaDistr level
960 Bool_t kFOUND = kFALSE;
961 for(Int_t ic = ncl[ir]; ic--;){
962 if(yres[ir][ic] - mean > 3. * syDis){
963 clst[ir][ic] = 0x0; continue;
964 }
965 nrow[nr]++; kFOUND = kTRUE;
966 }
967 // exit loop
968 if(kFOUND) nr++;
969 lr = ir; if(nr>=3) break;
29b87567 970 }
b1957d3c 971 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]);
972
973 // classify cluster rows
974 Int_t row = -1;
975 switch(nr){
976 case 1:
977 row = lr;
978 break;
979 case 2:
980 SetBit(kRowCross, kTRUE); // mark pad row crossing
981 if(nrow[0] > nrow[1]){ row = lr+1; lr = -1;}
982 else{
983 row = lr; lr = 1;
984 nrow[2] = nrow[1];
985 nrow[1] = nrow[0];
986 nrow[0] = nrow[2];
29b87567 987 }
b1957d3c 988 break;
989 case 3:
990 SetBit(kRowCross, kTRUE); // mark pad row crossing
991 break;
29b87567 992 }
b1957d3c 993 if(kPRINT) printf("\trow[%d] n[%d]\n\n", row, nrow[0]);
994 if(row<0) return kFALSE;
29b87567 995
b1957d3c 996 // Select and store clusters
997 // We should consider here :
998 // 1. How far is the chamber boundary
999 // 2. How big is the mean
3e778975 1000 Int_t n = 0;
b1957d3c 1001 for (Int_t ir = 0; ir < nr; ir++) {
1002 Int_t jr = row + ir*lr;
1003 if(kPRINT) printf("\tattach %d clusters for row %d\n", ncl[jr], jr);
1004 for (Int_t ic = 0; ic < ncl[jr]; ic++) {
1005 if(!(c = clst[jr][ic])) continue;
1006 Int_t it = c->GetPadTime();
1007 // TODO proper indexing of clusters !!
e3cf3d02 1008 fIndexes[it+kNtb*ir] = chamber->GetTB(it)->GetGlobalIndex(idxs[jr][ic]);
1009 fClusters[it+kNtb*ir] = c;
29b87567 1010
b1957d3c 1011 //printf("\tid[%2d] it[%d] idx[%d]\n", ic, it, fIndexes[it]);
1012
3e778975 1013 n++;
b1957d3c 1014 }
1015 }
1016
29b87567 1017 // number of minimum numbers of clusters expected for the tracklet
3e778975 1018 if (n < kClmin){
1019 //AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", n, kClmin));
e4f2f73d 1020 return kFALSE;
1021 }
3e778975 1022 SetN(n);
0906e73e 1023
e3cf3d02 1024 // Load calibration parameters for this tracklet
1025 Calibrate();
b1957d3c 1026
1027 // calculate dx for time bins in the drift region (calibration aware)
e3cf3d02 1028 Int_t irp = 0; Float_t x[2] = {0.,0.}; Int_t tb[2]={0,0};
b1957d3c 1029 for (Int_t it = t0; it < AliTRDtrackerV1::GetNTimeBins(); it++) {
1030 if(!fClusters[it]) continue;
1031 x[irp] = fClusters[it]->GetX();
1032 tb[irp] = it;
1033 irp++;
1034 if(irp==2) break;
e3cf3d02 1035 }
d86ed84c 1036 Int_t dtb = tb[1] - tb[0];
1037 fdX = dtb ? (x[0] - x[1]) / dtb : 0.15;
b1957d3c 1038
29b87567 1039 return kTRUE;
e4f2f73d 1040}
1041
03cef9b2 1042//____________________________________________________________
1043void AliTRDseedV1::Bootstrap(const AliTRDReconstructor *rec)
1044{
1045// Fill in all derived information. It has to be called after recovery from file or HLT.
1046// The primitive data are
1047// - list of clusters
1048// - detector (as the detector will be removed from clusters)
1049// - position of anode wire (fX0) - temporary
1050// - track reference position and direction
1051// - momentum of the track
1052// - time bin length [cm]
1053//
1054// A.Bercuci <A.Bercuci@gsi.de> Oct 30th 2008
1055//
1056 fReconstructor = rec;
1057 AliTRDgeometry g;
1058 AliTRDpadPlane *pp = g.GetPadPlane(fDet);
dd8059a8 1059 fPad[0] = pp->GetLengthIPad();
1060 fPad[1] = pp->GetWidthIPad();
1061 fPad[3] = TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle());
e3cf3d02 1062 //fSnp = fYref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);
1063 //fTgl = fZref[1];
3e778975 1064 Int_t n = 0, nshare = 0, nused = 0;
03cef9b2 1065 AliTRDcluster **cit = &fClusters[0];
8d2bec9e 1066 for(Int_t ic = kNclusters; ic--; cit++){
03cef9b2 1067 if(!(*cit)) return;
3e778975 1068 n++;
1069 if((*cit)->IsShared()) nshare++;
1070 if((*cit)->IsUsed()) nused++;
03cef9b2 1071 }
3e778975 1072 SetN(n); SetNUsed(nused); SetNShared(nshare);
e3cf3d02 1073 Fit();
03cef9b2 1074 CookLabels();
1075 GetProbability();
1076}
1077
1078
e4f2f73d 1079//____________________________________________________________________
d937ad7a 1080Bool_t AliTRDseedV1::Fit(Bool_t tilt, Int_t errors)
e4f2f73d 1081{
1082 //
1083 // Linear fit of the tracklet
1084 //
1085 // Parameters :
1086 //
1087 // Output :
1088 // True if successful
1089 //
1090 // Detailed description
1091 // 2. Check if tracklet crosses pad row boundary
1092 // 1. Calculate residuals in the y (r-phi) direction
1093 // 3. Do a Least Square Fit to the data
1094 //
1095
e3cf3d02 1096 if(!IsCalibrated()){
1097 AliWarning("Tracklet fit failed. Call Calibrate().");
1098 return kFALSE;
1099 }
1100
29b87567 1101 const Int_t kClmin = 8;
010d62b0 1102
9462866a 1103
1104 // cluster error parametrization parameters
010d62b0 1105 // 1. sy total charge
9462866a 1106 const Float_t sq0inv = 0.019962; // [1/q0]
1107 const Float_t sqb = 1.0281564; //[cm]
010d62b0 1108 // 2. sy for the PRF
1109 const Float_t scy[AliTRDgeometry::kNlayer][4] = {
d937ad7a 1110 {2.827e-02, 9.600e-04, 4.296e-01, 2.271e-02},
1111 {2.952e-02,-2.198e-04, 4.146e-01, 2.339e-02},
1112 {3.090e-02, 1.514e-03, 4.020e-01, 2.402e-02},
1113 {3.260e-02,-2.037e-03, 3.946e-01, 2.509e-02},
1114 {3.439e-02,-3.601e-04, 3.883e-01, 2.623e-02},
1115 {3.510e-02, 2.066e-03, 3.651e-01, 2.588e-02},
010d62b0 1116 };
1117 // 3. sy parallel to the track
d937ad7a 1118 const Float_t sy0 = 2.649e-02; // [cm]
1119 const Float_t sya = -8.864e-04; // [cm]
1120 const Float_t syb = -2.435e-01; // [cm]
1121
010d62b0 1122 // 4. sx parallel to the track
d937ad7a 1123 const Float_t sxgc = 5.427e-02;
1124 const Float_t sxgm = 7.783e-01;
1125 const Float_t sxgs = 2.743e-01;
1126 const Float_t sxe0 =-2.065e+00;
1127 const Float_t sxe1 =-2.978e-02;
1128
010d62b0 1129 // 5. sx perpendicular to the track
d937ad7a 1130// const Float_t sxd0 = 1.881e-02;
1131// const Float_t sxd1 =-4.101e-01;
1132// const Float_t sxd2 = 1.572e+00;
1133
2f7d6ac8 1134 // get track direction
1135 Double_t y0 = fYref[0];
1136 Double_t dydx = fYref[1];
1137 Double_t z0 = fZref[0];
1138 Double_t dzdx = fZref[1];
1139 Double_t yt, zt;
ae4e8b84 1140
e3cf3d02 1141 // calculation of tg^2(phi - a_L) and tg^2(a_L)
1142 Double_t tgg = (dydx-fExB)/(1.+dydx*fExB); tgg *= tgg;
1143 //Double_t exb2= fExB*fExB;
1144
b1957d3c 1145 //AliTRDtrackerV1::AliTRDLeastSquare fitterZ;
24d8660e 1146 TLinearFitter fitterY(1, "pol1");
29b87567 1147 // convertion factor from square to gauss distribution for sigma
b1957d3c 1148 //Double_t convert = 1./TMath::Sqrt(12.);
ae4e8b84 1149
29b87567 1150 // book cluster information
8d2bec9e 1151 Double_t qc[kNclusters], xc[kNclusters], yc[kNclusters], zc[kNclusters], sy[kNclusters];
e3cf3d02 1152
010d62b0 1153 Int_t ily = AliTRDgeometry::GetLayer(fDet);
dd8059a8 1154 Int_t n = 0;
9eb2d46c 1155 AliTRDcluster *c=0x0, **jc = &fClusters[0];
9eb2d46c 1156 for (Int_t ic=0; ic<kNtb; ic++, ++jc) {
b1957d3c 1157 //zRow[ic] = -1;
29b87567 1158 xc[ic] = -1.;
1159 yc[ic] = 999.;
1160 zc[ic] = 999.;
1161 sy[ic] = 0.;
b1957d3c 1162 //sz[ic] = 0.;
9eb2d46c 1163 if(!(c = (*jc))) continue;
29b87567 1164 if(!c->IsInChamber()) continue;
9462866a 1165
29b87567 1166 Float_t w = 1.;
1167 if(c->GetNPads()>4) w = .5;
1168 if(c->GetNPads()>5) w = .2;
010d62b0 1169
b1957d3c 1170 //zRow[fN] = c->GetPadRow();
dd8059a8 1171 qc[n] = TMath::Abs(c->GetQ());
d937ad7a 1172 // correct cluster position for PRF and v drift
e3cf3d02 1173 //Int_t jc = TMath::Max(fN-3, 0);
1174 //xc[fN] = c->GetXloc(fT0, fVD, &qc[jc], &xc[jc]/*, z0 - c->GetX()*dzdx*/);
1175 //Double_t s2 = fS2PRF + fDiffL*fDiffL*xc[fN]/(1.+2.*exb2)+tgg*xc[fN]*xc[fN]*exb2/12.;
dd8059a8 1176 //yc[fN] = c->GetYloc(s2, GetPadWidth(), xc[fN], fExB);
e3cf3d02 1177
1178 // uncalibrated cluster correction
1179 // TODO remove
d937ad7a 1180 Double_t x, y; GetClusterXY(c, x, y);
dd8059a8 1181 xc[n] = fX0 - x;
1182 yc[n] = y;
1183 zc[n] = c->GetZ();
2f7d6ac8 1184
1185 // extrapolated y value for the track
dd8059a8 1186 yt = y0 - xc[n]*dydx;
2f7d6ac8 1187 // extrapolated z value for the track
dd8059a8 1188 zt = z0 - xc[n]*dzdx;
2f7d6ac8 1189 // tilt correction
dd8059a8 1190 if(tilt) yc[n] -= GetTilt()*(zc[n] - zt);
2f7d6ac8 1191
010d62b0 1192 // ELABORATE CLUSTER ERROR
1193 // TODO to be moved to AliTRDcluster
010d62b0 1194 // basic y error (|| to track).
dd8059a8 1195 sy[n] = xc[n] < AliTRDgeometry::CamHght() ? 2. : sy0 + sya*TMath::Exp(1./(xc[n]+syb));
d937ad7a 1196 //printf("cluster[%d]\n\tsy[0] = %5.3e [um]\n", fN, sy[fN]*1.e4);
010d62b0 1197 // y error due to total charge
dd8059a8 1198 sy[n] += sqb*(1./qc[n] - sq0inv);
d937ad7a 1199 //printf("\tsy[1] = %5.3e [um]\n", sy[fN]*1.e4);
010d62b0 1200 // y error due to PRF
dd8059a8 1201 sy[n] += scy[ily][0]*TMath::Gaus(c->GetCenter(), scy[ily][1], scy[ily][2]) - scy[ily][3];
d937ad7a 1202 //printf("\tsy[2] = %5.3e [um]\n", sy[fN]*1.e4);
1203
dd8059a8 1204 sy[n] *= sy[n];
010d62b0 1205
1206 // ADD ERROR ON x
9462866a 1207 // error of drift length parallel to the track
dd8059a8 1208 Double_t sx = sxgc*TMath::Gaus(xc[n], sxgm, sxgs) + TMath::Exp(sxe0+sxe1*xc[n]); // [cm]
d937ad7a 1209 //printf("\tsx[0] = %5.3e [um]\n", sx*1.e4);
9462866a 1210 // error of drift length perpendicular to the track
1211 //sx += sxd0 + sxd1*d + sxd2*d*d;
d937ad7a 1212 sx *= sx; // square sx
d937ad7a 1213
9462866a 1214 // add error from ExB
dd8059a8 1215 if(errors>0) sy[n] += fExB*fExB*sx;
d937ad7a 1216 //printf("\tsy[3] = %5.3e [um^2]\n", sy[fN]*1.e8);
1217
1218 // global radial error due to misalignment/miscalibration
1219 Double_t sx0 = 0.; sx0 *= sx0;
1220 // add sx contribution to sy due to track angle
dd8059a8 1221 if(errors>1) sy[n] += tgg*(sx+sx0);
d937ad7a 1222 // TODO we should add tilt pad correction here
1223 //printf("\tsy[4] = %5.3e [um^2]\n", sy[fN]*1.e8);
dd8059a8 1224 c->SetSigmaY2(sy[n]);
d937ad7a 1225
dd8059a8 1226 sy[n] = TMath::Sqrt(sy[n]);
1227 fitterY.AddPoint(&xc[n], yc[n], sy[n]);
1228 n++;
29b87567 1229 }
47d5d320 1230 // to few clusters
dd8059a8 1231 if (n < kClmin) return kFALSE;
2f7d6ac8 1232
d937ad7a 1233 // fit XY
2f7d6ac8 1234 fitterY.Eval();
d937ad7a 1235 fYfit[0] = fitterY.GetParameter(0);
1236 fYfit[1] = -fitterY.GetParameter(1);
1237 // store covariance
1238 Double_t *p = fitterY.GetCovarianceMatrix();
1239 fCov[0] = p[0]; // variance of y0
1240 fCov[1] = p[1]; // covariance of y0, dydx
1241 fCov[2] = p[3]; // variance of dydx
b1957d3c 1242 // the ref radial position is set at the minimum of
1243 // the y variance of the tracklet
f29f13a6 1244 fX = -fCov[1]/fCov[2]; //fXref = fX0 - fXref;
1245 fS2Y = fCov[0] +2.*fX*fCov[1] + fX*fX*fCov[2];
b1957d3c 1246
1247 // fit XZ
1248 if(IsRowCross()){
1249 // TODO pad row cross position estimation !!!
1250 //AliInfo(Form("Padrow cross in detector %d", fDet));
dd8059a8 1251 fZfit[0] = .5*(zc[0]+zc[n-1]); fZfit[1] = 0.;
f29f13a6 1252 fS2Z = 0.02+1.55*fZref[1]; fS2Z *= fS2Z;
b1957d3c 1253 } else {
1254 fZfit[0] = zc[0]; fZfit[1] = 0.;
dd8059a8 1255 fS2Z = GetPadLength()*GetPadLength()/12.;
29b87567 1256 }
1257
29b87567 1258
b1957d3c 1259// // determine z offset of the fit
1260// Float_t zslope = 0.;
1261// Int_t nchanges = 0, nCross = 0;
1262// if(nz==2){ // tracklet is crossing pad row
1263// // Find the break time allowing one chage on pad-rows
1264// // with maximal number of accepted clusters
1265// Int_t padRef = zRow[0];
1266// for (Int_t ic=1; ic<fN; ic++) {
1267// if(zRow[ic] == padRef) continue;
1268//
1269// // debug
1270// if(zRow[ic-1] == zRow[ic]){
1271// printf("ERROR in pad row change!!!\n");
1272// }
1273//
1274// // evaluate parameters of the crossing point
1275// Float_t sx = (xc[ic-1] - xc[ic])*convert;
1276// fCross[0] = .5 * (xc[ic-1] + xc[ic]);
1277// fCross[2] = .5 * (zc[ic-1] + zc[ic]);
1278// fCross[3] = TMath::Max(dzdx * sx, .01);
1279// zslope = zc[ic-1] > zc[ic] ? 1. : -1.;
1280// padRef = zRow[ic];
1281// nCross = ic;
1282// nchanges++;
1283// }
1284// }
1285//
1286// // condition on nCross and reset nchanges TODO
1287//
1288// if(nchanges==1){
1289// if(dzdx * zslope < 0.){
1290// AliInfo("Tracklet-Track mismatch in dzdx. TODO.");
1291// }
1292//
1293//
1294// //zc[nc] = fitterZ.GetFunctionParameter(0);
1295// fCross[1] = fYfit[0] - fCross[0] * fYfit[1];
1296// fCross[0] = fX0 - fCross[0];
1297// }
29b87567 1298
29b87567 1299 return kTRUE;
e4f2f73d 1300}
1301
e4f2f73d 1302
f29f13a6 1303/*
e3cf3d02 1304//_____________________________________________________________________________
1305void AliTRDseedV1::FitMI()
1306{
1307//
1308// Fit the seed.
1309// Marian Ivanov's version
1310//
1311// linear fit on the y direction with respect to the reference direction.
1312// The residuals for each x (x = xc - x0) are deduced from:
1313// dy = y - yt (1)
1314// the tilting correction is written :
1315// y = yc + h*(zc-zt) (2)
1316// yt = y0+dy/dx*x (3)
1317// zt = z0+dz/dx*x (4)
1318// from (1),(2),(3) and (4)
1319// dy = yc - y0 - (dy/dx + h*dz/dx)*x + h*(zc-z0)
1320// the last term introduces the correction on y direction due to tilting pads. There are 2 ways to account for this:
1321// 1. use tilting correction for calculating the y
1322// 2. neglect tilting correction here and account for it in the error parametrization of the tracklet.
1323 const Float_t kRatio = 0.8;
1324 const Int_t kClmin = 5;
1325 const Float_t kmaxtan = 2;
1326
1327 if (TMath::Abs(fYref[1]) > kmaxtan){
1328 //printf("Exit: Abs(fYref[1]) = %3.3f, kmaxtan = %3.3f\n", TMath::Abs(fYref[1]), kmaxtan);
1329 return; // Track inclined too much
1330 }
1331
1332 Float_t sigmaexp = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction
dd8059a8 1333 Float_t ycrosscor = GetPadLength() * GetTilt() * 0.5; // Y correction for crossing
e3cf3d02 1334 Int_t fNChange = 0;
1335
1336 Double_t sumw;
1337 Double_t sumwx;
1338 Double_t sumwx2;
1339 Double_t sumwy;
1340 Double_t sumwxy;
1341 Double_t sumwz;
1342 Double_t sumwxz;
1343
1344 // Buffering: Leave it constant fot Performance issues
1345 Int_t zints[kNtb]; // Histograming of the z coordinate
1346 // Get 1 and second max probable coodinates in z
1347 Int_t zouts[2*kNtb];
1348 Float_t allowedz[kNtb]; // Allowed z for given time bin
1349 Float_t yres[kNtb]; // Residuals from reference
dd8059a8 1350 //Float_t anglecor = GetTilt() * fZref[1]; // Correction to the angle
e3cf3d02 1351
1352 Float_t pos[3*kNtb]; memset(pos, 0, 3*kNtb*sizeof(Float_t));
1353 Float_t *fX = &pos[0], *fY = &pos[kNtb], *fZ = &pos[2*kNtb];
1354
1355 Int_t fN = 0; AliTRDcluster *c = 0x0;
1356 fN2 = 0;
1357 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1358 yres[i] = 10000.0;
1359 if (!(c = fClusters[i])) continue;
1360 if(!c->IsInChamber()) continue;
1361 // Residual y
dd8059a8 1362 //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] + GetTilt()*(fZ[i] - fZref[0]);
e3cf3d02 1363 fX[i] = fX0 - c->GetX();
1364 fY[i] = c->GetY();
1365 fZ[i] = c->GetZ();
dd8059a8 1366 yres[i] = fY[i] - GetTilt()*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
e3cf3d02 1367 zints[fN] = Int_t(fZ[i]);
1368 fN++;
1369 }
1370
1371 if (fN < kClmin){
1372 //printf("Exit fN < kClmin: fN = %d\n", fN);
1373 return;
1374 }
1375 Int_t nz = AliTRDtrackerV1::Freq(fN, zints, zouts, kFALSE);
1376 Float_t fZProb = zouts[0];
1377 if (nz <= 1) zouts[3] = 0;
1378 if (zouts[1] + zouts[3] < kClmin) {
1379 //printf("Exit zouts[1] = %d, zouts[3] = %d\n",zouts[1],zouts[3]);
1380 return;
1381 }
1382
1383 // Z distance bigger than pad - length
1384 if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) zouts[3] = 0;
1385
1386 Int_t breaktime = -1;
1387 Bool_t mbefore = kFALSE;
1388 Int_t cumul[kNtb][2];
1389 Int_t counts[2] = { 0, 0 };
1390
1391 if (zouts[3] >= 3) {
1392
1393 //
1394 // Find the break time allowing one chage on pad-rows
1395 // with maximal number of accepted clusters
1396 //
1397 fNChange = 1;
1398 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1399 cumul[i][0] = counts[0];
1400 cumul[i][1] = counts[1];
1401 if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++;
1402 if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++;
1403 }
1404 Int_t maxcount = 0;
1405 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1406 Int_t after = cumul[AliTRDtrackerV1::GetNTimeBins()][0] - cumul[i][0];
1407 Int_t before = cumul[i][1];
1408 if (after + before > maxcount) {
1409 maxcount = after + before;
1410 breaktime = i;
1411 mbefore = kFALSE;
1412 }
1413 after = cumul[AliTRDtrackerV1::GetNTimeBins()-1][1] - cumul[i][1];
1414 before = cumul[i][0];
1415 if (after + before > maxcount) {
1416 maxcount = after + before;
1417 breaktime = i;
1418 mbefore = kTRUE;
1419 }
1420 }
1421 breaktime -= 1;
1422 }
1423
1424 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1425 if (i > breaktime) allowedz[i] = mbefore ? zouts[2] : zouts[0];
1426 if (i <= breaktime) allowedz[i] = (!mbefore) ? zouts[2] : zouts[0];
1427 }
1428
1429 if (((allowedz[0] > allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] < 0)) ||
1430 ((allowedz[0] < allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] > 0))) {
1431 //
1432 // Tracklet z-direction not in correspondance with track z direction
1433 //
1434 fNChange = 0;
1435 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1436 allowedz[i] = zouts[0]; // Only longest taken
1437 }
1438 }
1439
1440 if (fNChange > 0) {
1441 //
1442 // Cross pad -row tracklet - take the step change into account
1443 //
1444 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1445 if (!fClusters[i]) continue;
1446 if(!fClusters[i]->IsInChamber()) continue;
1447 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
1448 // Residual y
dd8059a8 1449 //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] + GetTilt()*(fZ[i] - fZref[0]);
1450 yres[i] = fY[i] - GetTilt()*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
f29f13a6 1451// if (TMath::Abs(fZ[i] - fZProb) > 2) {
dd8059a8 1452// if (fZ[i] > fZProb) yres[i] += GetTilt() * GetPadLength();
1453// if (fZ[i] < fZProb) yres[i] -= GetTilt() * GetPadLength();
f29f13a6 1454 }
e3cf3d02 1455 }
1456 }
1457
1458 Double_t yres2[kNtb];
1459 Double_t mean;
1460 Double_t sigma;
1461 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1462 if (!fClusters[i]) continue;
1463 if(!fClusters[i]->IsInChamber()) continue;
1464 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
1465 yres2[fN2] = yres[i];
1466 fN2++;
1467 }
1468 if (fN2 < kClmin) {
1469 //printf("Exit fN2 < kClmin: fN2 = %d\n", fN2);
1470 fN2 = 0;
1471 return;
1472 }
1473 AliMathBase::EvaluateUni(fN2,yres2,mean,sigma, Int_t(fN2*kRatio-2.));
1474 if (sigma < sigmaexp * 0.8) {
1475 sigma = sigmaexp;
1476 }
1477 //Float_t fSigmaY = sigma;
1478
1479 // Reset sums
1480 sumw = 0;
1481 sumwx = 0;
1482 sumwx2 = 0;
1483 sumwy = 0;
1484 sumwxy = 0;
1485 sumwz = 0;
1486 sumwxz = 0;
1487
1488 fN2 = 0;
1489 Float_t fMeanz = 0;
1490 Float_t fMPads = 0;
1491 fUsable = 0;
1492 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1493 if (!fClusters[i]) continue;
1494 if (!fClusters[i]->IsInChamber()) continue;
1495 if (TMath::Abs(fZ[i] - allowedz[i]) > 2){fClusters[i] = 0x0; continue;}
1496 if (TMath::Abs(yres[i] - mean) > 4.0 * sigma){fClusters[i] = 0x0; continue;}
1497 SETBIT(fUsable,i);
1498 fN2++;
1499 fMPads += fClusters[i]->GetNPads();
1500 Float_t weight = 1.0;
1501 if (fClusters[i]->GetNPads() > 4) weight = 0.5;
1502 if (fClusters[i]->GetNPads() > 5) weight = 0.2;
1503
1504
1505 Double_t x = fX[i];
1506 //printf("x = %7.3f dy = %7.3f fit %7.3f\n", x, yres[i], fY[i]-yres[i]);
1507
1508 sumw += weight;
1509 sumwx += x * weight;
1510 sumwx2 += x*x * weight;
1511 sumwy += weight * yres[i];
1512 sumwxy += weight * (yres[i]) * x;
1513 sumwz += weight * fZ[i];
1514 sumwxz += weight * fZ[i] * x;
1515
1516 }
1517
1518 if (fN2 < kClmin){
1519 //printf("Exit fN2 < kClmin(2): fN2 = %d\n",fN2);
1520 fN2 = 0;
1521 return;
1522 }
1523 fMeanz = sumwz / sumw;
1524 Float_t correction = 0;
1525 if (fNChange > 0) {
1526 // Tracklet on boundary
1527 if (fMeanz < fZProb) correction = ycrosscor;
1528 if (fMeanz > fZProb) correction = -ycrosscor;
1529 }
1530
1531 Double_t det = sumw * sumwx2 - sumwx * sumwx;
1532 fYfit[0] = (sumwx2 * sumwy - sumwx * sumwxy) / det;
1533 fYfit[1] = (sumw * sumwxy - sumwx * sumwy) / det;
1534
1535 fS2Y = 0;
1536 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1537 if (!TESTBIT(fUsable,i)) continue;
1538 Float_t delta = yres[i] - fYfit[0] - fYfit[1] * fX[i];
1539 fS2Y += delta*delta;
1540 }
1541 fS2Y = TMath::Sqrt(fS2Y / Float_t(fN2-2));
1542 // TEMPORARY UNTIL covariance properly calculated
1543 fS2Y = TMath::Max(fS2Y, Float_t(.1));
1544
1545 fZfit[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det;
1546 fZfit[1] = (sumw * sumwxz - sumwx * sumwz) / det;
1547// fYfitR[0] += fYref[0] + correction;
1548// fYfitR[1] += fYref[1];
1549// fYfit[0] = fYfitR[0];
1550 fYfit[1] = -fYfit[1];
1551
1552 UpdateUsed();
f29f13a6 1553}*/
e3cf3d02 1554
e4f2f73d 1555//___________________________________________________________________
203967fc 1556void AliTRDseedV1::Print(Option_t *o) const
e4f2f73d 1557{
1558 //
1559 // Printing the seedstatus
1560 //
1561
eb2b4f91 1562 AliInfo(Form("Det[%3d] X0[%7.2f] Pad[L[%5.2f] W[%5.2f] Tilt[%+6.2f]]", fDet, fX0, GetPadLength(), GetPadWidth(), GetTilt()));
dd8059a8 1563 AliInfo(Form("N[%2d] Nused[%2d] Nshared[%2d] [%d]", GetN(), GetNUsed(), GetNShared(), fN));
1564
1565 Double_t cov[3], x=GetX();
1566 GetCovAt(x, cov);
1567 AliInfo(" | x[cm] | y[cm] | z[cm] | dydx | dzdx |");
1568 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]));
1569 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[2]), fZref[0]-fX*fYref[1], TMath::Sqrt(fRefCov[2]), fYref[1], fZref[1]))
203967fc 1570
1571
1572 if(strcmp(o, "a")!=0) return;
1573
4dc4dc2e 1574 AliTRDcluster* const* jc = &fClusters[0];
8d2bec9e 1575 for(int ic=0; ic<kNclusters; ic++, jc++) {
4dc4dc2e 1576 if(!(*jc)) continue;
203967fc 1577 (*jc)->Print(o);
4dc4dc2e 1578 }
e4f2f73d 1579}
47d5d320 1580
203967fc 1581
1582//___________________________________________________________________
1583Bool_t AliTRDseedV1::IsEqual(const TObject *o) const
1584{
1585 // Checks if current instance of the class has the same essential members
1586 // as the given one
1587
1588 if(!o) return kFALSE;
1589 const AliTRDseedV1 *inTracklet = dynamic_cast<const AliTRDseedV1*>(o);
1590 if(!inTracklet) return kFALSE;
1591
1592 for (Int_t i = 0; i < 2; i++){
e3cf3d02 1593 if ( fYref[i] != inTracklet->fYref[i] ) return kFALSE;
1594 if ( fZref[i] != inTracklet->fZref[i] ) return kFALSE;
203967fc 1595 }
1596
e3cf3d02 1597 if ( fS2Y != inTracklet->fS2Y ) return kFALSE;
dd8059a8 1598 if ( GetTilt() != inTracklet->GetTilt() ) return kFALSE;
1599 if ( GetPadLength() != inTracklet->GetPadLength() ) return kFALSE;
203967fc 1600
8d2bec9e 1601 for (Int_t i = 0; i < kNclusters; i++){
e3cf3d02 1602// if ( fX[i] != inTracklet->GetX(i) ) return kFALSE;
1603// if ( fY[i] != inTracklet->GetY(i) ) return kFALSE;
1604// if ( fZ[i] != inTracklet->GetZ(i) ) return kFALSE;
1605 if ( fIndexes[i] != inTracklet->fIndexes[i] ) return kFALSE;
203967fc 1606 }
f29f13a6 1607// if ( fUsable != inTracklet->fUsable ) return kFALSE;
203967fc 1608
1609 for (Int_t i=0; i < 2; i++){
e3cf3d02 1610 if ( fYfit[i] != inTracklet->fYfit[i] ) return kFALSE;
1611 if ( fZfit[i] != inTracklet->fZfit[i] ) return kFALSE;
1612 if ( fLabels[i] != inTracklet->fLabels[i] ) return kFALSE;
203967fc 1613 }
1614
e3cf3d02 1615/* if ( fMeanz != inTracklet->GetMeanz() ) return kFALSE;
1616 if ( fZProb != inTracklet->GetZProb() ) return kFALSE;*/
3e778975 1617 if ( fN != inTracklet->fN ) return kFALSE;
1618 //if ( fNUsed != inTracklet->fNUsed ) return kFALSE;
e3cf3d02 1619 //if ( fFreq != inTracklet->GetFreq() ) return kFALSE;
1620 //if ( fNChange != inTracklet->GetNChange() ) return kFALSE;
203967fc 1621
e3cf3d02 1622 if ( fC != inTracklet->fC ) return kFALSE;
1623 //if ( fCC != inTracklet->GetCC() ) return kFALSE;
1624 if ( fChi2 != inTracklet->fChi2 ) return kFALSE;
203967fc 1625 // if ( fChi2Z != inTracklet->GetChi2Z() ) return kFALSE;
1626
e3cf3d02 1627 if ( fDet != inTracklet->fDet ) return kFALSE;
1628 if ( fMom != inTracklet->fMom ) return kFALSE;
1629 if ( fdX != inTracklet->fdX ) return kFALSE;
203967fc 1630
8d2bec9e 1631 for (Int_t iCluster = 0; iCluster < kNclusters; iCluster++){
203967fc 1632 AliTRDcluster *curCluster = fClusters[iCluster];
e3cf3d02 1633 AliTRDcluster *inCluster = inTracklet->fClusters[iCluster];
203967fc 1634 if (curCluster && inCluster){
1635 if (! curCluster->IsEqual(inCluster) ) {
1636 curCluster->Print();
1637 inCluster->Print();
1638 return kFALSE;
1639 }
1640 } else {
1641 // if one cluster exists, and corresponding
1642 // in other tracklet doesn't - return kFALSE
1643 if(curCluster || inCluster) return kFALSE;
1644 }
1645 }
1646 return kTRUE;
1647}