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