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