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