Using SetMaxDCAToVertexXY instead of the removed SetDCAToVertex
[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////////////////////////////////////////////////////////////////////////////
19// //
20// The TRD track seed //
21// //
22// Authors: //
23// Alex Bercuci <A.Bercuci@gsi.de> //
24// Markus Fasel <M.Fasel@gsi.de> //
25// //
26////////////////////////////////////////////////////////////////////////////
27
28#include "TMath.h"
29#include "TLinearFitter.h"
eb38ed55 30#include "TClonesArray.h" // tmp
31#include <TTreeStream.h>
e4f2f73d 32
33#include "AliLog.h"
34#include "AliMathBase.h"
d937ad7a 35#include "AliCDBManager.h"
36#include "AliTracker.h"
e4f2f73d 37
03cef9b2 38#include "AliTRDpadPlane.h"
e4f2f73d 39#include "AliTRDcluster.h"
f3d3af1b 40#include "AliTRDseedV1.h"
41#include "AliTRDtrackV1.h"
e4f2f73d 42#include "AliTRDcalibDB.h"
eb38ed55 43#include "AliTRDchamberTimeBin.h"
44#include "AliTRDtrackingChamber.h"
45#include "AliTRDtrackerV1.h"
46#include "AliTRDReconstructor.h"
e4f2f73d 47#include "AliTRDrecoParam.h"
a076fc2f 48#include "AliTRDCommonParam.h"
d937ad7a 49
0906e73e 50#include "Cal/AliTRDCalPID.h"
d937ad7a 51#include "Cal/AliTRDCalROC.h"
52#include "Cal/AliTRDCalDet.h"
e4f2f73d 53
e4f2f73d 54ClassImp(AliTRDseedV1)
55
56//____________________________________________________________________
ae4e8b84 57AliTRDseedV1::AliTRDseedV1(Int_t det)
e4f2f73d 58 :AliTRDseed()
3a039a31 59 ,fReconstructor(0x0)
ae4e8b84 60 ,fClusterIter(0x0)
61 ,fClusterIdx(0)
62 ,fDet(det)
0906e73e 63 ,fMom(0.)
bcb6fb78 64 ,fSnp(0.)
65 ,fTgl(0.)
66 ,fdX(0.)
6e4d4425 67 ,fXref(0.)
d937ad7a 68 ,fExB(0.)
e4f2f73d 69{
70 //
71 // Constructor
72 //
29b87567 73 for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = 0.;
74 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = -1.;
6e4d4425 75 fRefCov[0] = 1.; fRefCov[1] = 0.; fRefCov[2] = 1.;
d937ad7a 76 // covariance matrix [diagonal]
77 // default sy = 200um and sz = 2.3 cm
78 fCov[0] = 4.e-4; fCov[1] = 0.; fCov[2] = 5.3;
e4f2f73d 79}
80
81//____________________________________________________________________
0906e73e 82AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
e4f2f73d 83 :AliTRDseed((AliTRDseed&)ref)
43d6ad34 84 ,fReconstructor(ref.fReconstructor)
ae4e8b84 85 ,fClusterIter(0x0)
86 ,fClusterIdx(0)
87 ,fDet(ref.fDet)
0906e73e 88 ,fMom(ref.fMom)
bcb6fb78 89 ,fSnp(ref.fSnp)
90 ,fTgl(ref.fTgl)
91 ,fdX(ref.fdX)
6e4d4425 92 ,fXref(ref.fXref)
d937ad7a 93 ,fExB(ref.fExB)
e4f2f73d 94{
95 //
96 // Copy Constructor performing a deep copy
97 //
98
29b87567 99 //printf("AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &)\n");
100 SetBit(kOwner, kFALSE);
101 for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = ref.fdEdx[islice];
102 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = ref.fProb[ispec];
6e4d4425 103 memcpy(fRefCov, ref.fRefCov, 3*sizeof(Double_t));
d937ad7a 104 memcpy(fCov, ref.fCov, 3*sizeof(Double_t));
fbb2ea06 105}
d9950a5a 106
0906e73e 107
e4f2f73d 108//____________________________________________________________________
109AliTRDseedV1& AliTRDseedV1::operator=(const AliTRDseedV1 &ref)
110{
111 //
112 // Assignment Operator using the copy function
113 //
114
29b87567 115 if(this != &ref){
116 ref.Copy(*this);
117 }
221ab7e0 118 SetBit(kOwner, kFALSE);
119
29b87567 120 return *this;
e4f2f73d 121
122}
123
124//____________________________________________________________________
125AliTRDseedV1::~AliTRDseedV1()
126{
127 //
128 // Destructor. The RecoParam object belongs to the underlying tracker.
129 //
130
29b87567 131 //printf("I-AliTRDseedV1::~AliTRDseedV1() : Owner[%s]\n", IsOwner()?"YES":"NO");
e4f2f73d 132
29b87567 133 if(IsOwner())
134 for(int itb=0; itb<knTimebins; itb++){
135 if(!fClusters[itb]) continue;
136 //AliInfo(Form("deleting c %p @ %d", fClusters[itb], itb));
137 delete fClusters[itb];
138 fClusters[itb] = 0x0;
139 }
e4f2f73d 140}
141
142//____________________________________________________________________
143void AliTRDseedV1::Copy(TObject &ref) const
144{
145 //
146 // Copy function
147 //
148
29b87567 149 //AliInfo("");
150 AliTRDseedV1 &target = (AliTRDseedV1 &)ref;
151
ae4e8b84 152 target.fClusterIter = 0x0;
153 target.fClusterIdx = 0;
154 target.fDet = fDet;
29b87567 155 target.fMom = fMom;
156 target.fSnp = fSnp;
157 target.fTgl = fTgl;
158 target.fdX = fdX;
6e4d4425 159 target.fXref = fXref;
d937ad7a 160 target.fExB = fExB;
29b87567 161 target.fReconstructor = fReconstructor;
162
163 for(int islice=0; islice < knSlices; islice++) target.fdEdx[islice] = fdEdx[islice];
164 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) target.fProb[ispec] = fProb[ispec];
6e4d4425 165 memcpy(target.fRefCov, fRefCov, 3*sizeof(Double_t));
d937ad7a 166 memcpy(target.fCov, fCov, 3*sizeof(Double_t));
29b87567 167
168 AliTRDseed::Copy(target);
e4f2f73d 169}
170
0906e73e 171
172//____________________________________________________________
f3d3af1b 173Bool_t AliTRDseedV1::Init(AliTRDtrackV1 *track)
0906e73e 174{
175// Initialize this tracklet using the track information
176//
177// Parameters:
178// track - the TRD track used to initialize the tracklet
179//
180// Detailed description
181// The function sets the starting point and direction of the
182// tracklet according to the information from the TRD track.
183//
184// Caution
185// The TRD track has to be propagated to the beginning of the
186// chamber where the tracklet will be constructed
187//
188
29b87567 189 Double_t y, z;
190 if(!track->GetProlongation(fX0, y, z)) return kFALSE;
b1957d3c 191 UpDate(track);
29b87567 192 return kTRUE;
0906e73e 193}
194
bcb6fb78 195
196//____________________________________________________________________
b1957d3c 197void AliTRDseedV1::UpDate(const AliTRDtrackV1 *trk)
198{
199 // update tracklet reference position from the TRD track
200 // Funny name to avoid the clash with the function AliTRDseed::Update() (has to be made obselete)
201
202 fSnp = trk->GetSnp();
203 fTgl = trk->GetTgl();
204 fMom = trk->GetP();
205 fYref[1] = fSnp/(1. - fSnp*fSnp);
206 fZref[1] = fTgl;
207 SetCovRef(trk->GetCovariance());
208
209 Double_t dx = trk->GetX() - fX0;
210 fYref[0] = trk->GetY() - dx*fYref[1];
211 fZref[0] = trk->GetZ() - dx*fZref[1];
212}
213
214//____________________________________________________________________
bcb6fb78 215void AliTRDseedV1::CookdEdx(Int_t nslices)
216{
217// Calculates average dE/dx for all slices and store them in the internal array fdEdx.
218//
219// Parameters:
220// nslices : number of slices for which dE/dx should be calculated
221// Output:
222// store results in the internal array fdEdx. This can be accessed with the method
223// AliTRDseedV1::GetdEdx()
224//
225// Detailed description
226// Calculates average dE/dx for all slices. Depending on the PID methode
227// the number of slices can be 3 (LQ) or 8(NN).
3ee48d6e 228// The calculation of dQ/dl are done using the tracklet fit results (see AliTRDseedV1::GetdQdl(Int_t))
bcb6fb78 229//
230// The following effects are included in the calculation:
231// 1. calibration values for t0 and vdrift (using x coordinate to calculate slice)
232// 2. cluster sharing (optional see AliTRDrecoParam::SetClusterSharing())
233// 3. cluster size
234//
235
29b87567 236 Int_t nclusters[knSlices];
237 for(int i=0; i<knSlices; i++){
238 fdEdx[i] = 0.;
239 nclusters[i] = 0;
240 }
e73abf77 241 const Double_t kDriftLength = (.5 * AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
29b87567 242
3ee48d6e 243 AliTRDcluster *c = 0x0;
29b87567 244 for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++){
8e709c82 245 if(!(c = fClusters[ic]) && !(c = fClusters[ic+kNtb])) continue;
e73abf77 246 Float_t dx = TMath::Abs(fX0 - c->GetX());
29b87567 247
248 // Filter clusters for dE/dx calculation
249
250 // 1.consider calibration effects for slice determination
e73abf77 251 Int_t slice;
252 if(dx<kDriftLength){ // TODO should be replaced by c->IsInChamber()
253 slice = Int_t(dx * nslices / kDriftLength);
254 } else slice = c->GetX() < fX0 ? nslices-1 : 0;
255
256
29b87567 257 // 2. take sharing into account
3ee48d6e 258 Float_t w = c->IsShared() ? .5 : 1.;
29b87567 259
260 // 3. take into account large clusters TODO
261 //w *= c->GetNPads() > 3 ? .8 : 1.;
262
263 //CHECK !!!
264 fdEdx[slice] += w * GetdQdl(ic); //fdQdl[ic];
265 nclusters[slice]++;
266 } // End of loop over clusters
267
cd40b287 268 //if(fReconstructor->GetPIDMethod() == AliTRDReconstructor::kLQPID){
0d83b3a5 269 if(nslices == AliTRDpidUtil::kLQslices){
29b87567 270 // calculate mean charge per slice (only LQ PID)
271 for(int is=0; is<nslices; is++){
272 if(nclusters[is]) fdEdx[is] /= nclusters[is];
273 }
274 }
bcb6fb78 275}
276
d937ad7a 277//____________________________________________________________________
278void AliTRDseedV1::GetClusterXY(const AliTRDcluster *c, Double_t &x, Double_t &y)
279{
280// Return corrected position of the cluster taking into
281// account variation of the drift velocity with drift length.
282
283
284 // drift velocity correction TODO to be moved to the clusterizer
285 const Float_t cx[] = {
286 -9.6280e-02, 1.3091e-01,-1.7415e-02,-9.9221e-02,-1.2040e-01,-9.5493e-02,
287 -5.0041e-02,-1.6726e-02, 3.5756e-03, 1.8611e-02, 2.6378e-02, 3.3823e-02,
288 3.4811e-02, 3.5282e-02, 3.5386e-02, 3.6047e-02, 3.5201e-02, 3.4384e-02,
289 3.2864e-02, 3.1932e-02, 3.2051e-02, 2.2539e-02,-2.5154e-02,-1.2050e-01,
290 -1.2050e-01
291 };
292
293 // PRF correction TODO to be replaced by the gaussian
294 // approximation with full error parametrization and // moved to the clusterizer
295 const Float_t cy[AliTRDgeometry::kNlayer][3] = {
296 { 4.014e-04, 8.605e-03, -6.880e+00},
297 {-3.061e-04, 9.663e-03, -6.789e+00},
298 { 1.124e-03, 1.105e-02, -6.825e+00},
299 {-1.527e-03, 1.231e-02, -6.777e+00},
300 { 2.150e-03, 1.387e-02, -6.783e+00},
301 {-1.296e-03, 1.486e-02, -6.825e+00}
302 };
303
304 Int_t ily = AliTRDgeometry::GetLayer(c->GetDetector());
305 x = c->GetX() - cx[c->GetLocalTimeBin()];
306 y = c->GetY() + cy[ily][0] + cy[ily][1] * TMath::Sin(cy[ily][2] * c->GetCenter());
307 return;
308}
b83573da 309
bcb6fb78 310//____________________________________________________________________
311Float_t AliTRDseedV1::GetdQdl(Int_t ic) const
312{
3ee48d6e 313// Using the linear approximation of the track inside one TRD chamber (TRD tracklet)
314// the charge per unit length can be written as:
315// BEGIN_LATEX
316// #frac{dq}{dl} = #frac{q_{c}}{dx * #sqrt{1 + #(){#frac{dy}{dx}}^{2}_{fit} + #(){#frac{dy}{dx}}^{2}_{ref}}}
317// END_LATEX
318// where qc is the total charge collected in the current time bin and dx is the length
319// of the time bin. For the moment (Jan 20 2009) only pad row cross corrections are
320// considered for the charge but none are applied for drift velocity variations along
321// the drift region or assymetry of the TRF
322//
323// Author : Alex Bercuci <A.Bercuci@gsi.de>
324//
325 Float_t dq = 0.;
326 if(fClusters[ic]) dq += TMath::Abs(fClusters[ic]->GetQ());
8e709c82 327 if(fClusters[ic+kNtb]) dq += TMath::Abs(fClusters[ic+kNtb]->GetQ());
328 if(dq<1.e-3 || fdX < 1.e-3) return 0.;
3ee48d6e 329
330 return dq/fdX/TMath::Sqrt(1. + fYfit[1]*fYfit[1] + fZref[1]*fZref[1]);
bcb6fb78 331}
332
0906e73e 333//____________________________________________________________________
334Double_t* AliTRDseedV1::GetProbability()
335{
336// Fill probability array for tracklet from the DB.
337//
338// Parameters
339//
340// Output
341// returns pointer to the probability array and 0x0 if missing DB access
342//
343// Detailed description
344
29b87567 345
346 // retrive calibration db
0906e73e 347 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
348 if (!calibration) {
349 AliError("No access to calibration data");
350 return 0x0;
351 }
352
3a039a31 353 if (!fReconstructor) {
354 AliError("Reconstructor not set.");
4ba1d6ae 355 return 0x0;
356 }
357
0906e73e 358 // Retrieve the CDB container class with the parametric detector response
3a039a31 359 const AliTRDCalPID *pd = calibration->GetPIDObject(fReconstructor->GetPIDMethod());
0906e73e 360 if (!pd) {
361 AliError("No access to AliTRDCalPID object");
362 return 0x0;
363 }
29b87567 364 //AliInfo(Form("Method[%d] : %s", fReconstructor->GetRecoParam() ->GetPIDMethod(), pd->IsA()->GetName()));
10f75631 365
29b87567 366 // calculate tracklet length TO DO
0906e73e 367 Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
368 /// TMath::Sqrt((1.0 - fSnp[iPlane]*fSnp[iPlane]) / (1.0 + fTgl[iPlane]*fTgl[iPlane]));
369
370 //calculate dE/dx
3a039a31 371 CookdEdx(fReconstructor->GetNdEdxSlices());
0906e73e 372
373 // Sets the a priori probabilities
374 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
ae4e8b84 375 fProb[ispec] = pd->GetProbability(ispec, fMom, &fdEdx[0], length, GetPlane());
0906e73e 376 }
377
29b87567 378 return &fProb[0];
0906e73e 379}
380
381//____________________________________________________________________
e4f2f73d 382Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
383{
384 //
385 // Returns a quality measurement of the current seed
386 //
387
29b87567 388 Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
389 return
390 .5 * TMath::Abs(18.0 - fN2)
391 + 10.* TMath::Abs(fYfit[1] - fYref[1])
392 + 5. * TMath::Abs(fYfit[0] - fYref[0] + zcorr)
393 + 2. * TMath::Abs(fMeanz - fZref[0]) / fPadLength;
e4f2f73d 394}
395
396//____________________________________________________________________
d937ad7a 397void AliTRDseedV1::GetCovAt(Double_t x, Double_t *cov) const
0906e73e 398{
d937ad7a 399// Computes covariance in the y-z plane at radial point x (in tracking coordinates)
400// and returns the results in the preallocated array cov[3] as :
401// cov[0] = Var(y)
402// cov[1] = Cov(yz)
403// cov[2] = Var(z)
404//
405// Details
406//
407// For the linear transformation
408// BEGIN_LATEX
409// Y = T_{x} X^{T}
410// END_LATEX
411// The error propagation has the general form
412// BEGIN_LATEX
413// C_{Y} = T_{x} C_{X} T_{x}^{T}
414// END_LATEX
415// We apply this formula 2 times. First to calculate the covariance of the tracklet
416// at point x we consider:
417// BEGIN_LATEX
418// T_{x} = (1 x); X=(y0 dy/dx); C_{X}=#(){#splitline{Var(y0) Cov(y0, dy/dx)}{Cov(y0, dy/dx) Var(dy/dx)}}
419// END_LATEX
420// and secondly to take into account the tilt angle
421// BEGIN_LATEX
422// T_{#alpha} = #(){#splitline{cos(#alpha) __ sin(#alpha)}{-sin(#alpha) __ cos(#alpha)}}; X=(y z); C_{X}=#(){#splitline{Var(y) 0}{0 Var(z)}}
423// END_LATEX
424//
425// using simple trigonometrics one can write for this last case
426// BEGIN_LATEX
427// 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})}}
428// END_LATEX
429// which can be aproximated for small alphas (2 deg) with
430// BEGIN_LATEX
431// 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}}}
432// END_LATEX
433//
434// before applying the tilt rotation we also apply systematic uncertainties to the tracklet
435// position which can be tunned from outside via the AliTRDrecoParam::SetSysCovMatrix(). They might
436// account for extra misalignment/miscalibration uncertainties.
437//
438// Author :
439// Alex Bercuci <A.Bercuci@gsi.de>
440// Date : Jan 8th 2009
441//
b1957d3c 442
443
d937ad7a 444 Double_t xr = fX0-x;
445 Double_t sy2 = fCov[0] +2.*xr*fCov[1] + xr*xr*fCov[2];
446 Double_t sz2 = fPadLength*fPadLength/12.;
0906e73e 447
d937ad7a 448 // insert systematic uncertainties
449 Double_t sys[15];
450 fReconstructor->GetRecoParam()->GetSysCovMatrix(sys);
451 sy2 += sys[0];
452 sz2 += sys[1];
453
454 // rotate covariance matrix
455 Double_t t2 = fTilt*fTilt;
456 Double_t correction = 1./(1. + t2);
457 cov[0] = (sy2+t2*sz2)*correction;
458 cov[1] = fTilt*(sz2 - sy2)*correction;
459 cov[2] = (t2*sy2+sz2)*correction;
460}
eb38ed55 461
0906e73e 462
d937ad7a 463//____________________________________________________________________
464void AliTRDseedV1::SetExB()
465{
466// Retrive the tg(a_L) from OCDB. The following information are used
467// - detector index
468// - column and row position of first attached cluster.
469//
470// If no clusters are attached to the tracklet a random central chamber
471// position (c=70, r=7) will be used to retrieve the drift velocity.
472//
473// Author :
474// Alex Bercuci <A.Bercuci@gsi.de>
475// Date : Jan 8th 2009
476//
eb38ed55 477
d937ad7a 478 AliCDBManager *cdb = AliCDBManager::Instance();
479 if(cdb->GetRun() < 0){
480 AliError("OCDB manager not properly initialized");
481 return;
482 }
0906e73e 483
d937ad7a 484 AliTRDcalibDB *fCalib = AliTRDcalibDB::Instance();
485 AliTRDCalROC *fCalVdriftROC = fCalib->GetVdriftROC(fDet);
486 const AliTRDCalDet *fCalVdriftDet = fCalib->GetVdriftDet();
487
488 Int_t col = 70, row = 7;
489 AliTRDcluster **c = &fClusters[0];
490 if(fN){
491 Int_t ic = 0;
492 while (ic<AliTRDseed::knTimebins && !(*c)){ic++; c++;}
493 if(*c){
494 col = (*c)->GetPadCol();
495 row = (*c)->GetPadRow();
496 }
497 }
3a039a31 498
d937ad7a 499 Double_t vd = fCalVdriftDet->GetValue(fDet) * fCalVdriftROC->GetValue(col, row);
a076fc2f 500 fExB = AliTRDCommonParam::Instance()->GetOmegaTau(vd);
0906e73e 501}
502
0906e73e 503//____________________________________________________________________
29b87567 504void AliTRDseedV1::SetOwner()
0906e73e 505{
29b87567 506 //AliInfo(Form("own [%s] fOwner[%s]", own?"YES":"NO", fOwner?"YES":"NO"));
507
508 if(TestBit(kOwner)) return;
509 for(int ic=0; ic<knTimebins; ic++){
510 if(!fClusters[ic]) continue;
511 fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
512 }
513 SetBit(kOwner);
0906e73e 514}
515
516//____________________________________________________________________
eb38ed55 517Bool_t AliTRDseedV1::AttachClustersIter(AliTRDtrackingChamber *chamber, Float_t quality, Bool_t kZcorr, AliTRDcluster *c)
e4f2f73d 518{
519 //
520 // Iterative process to register clusters to the seed.
521 // In iteration 0 we try only one pad-row and if quality not
522 // sufficient we try 2 pad-rows (about 5% of tracks cross 2 pad-rows)
523 //
29b87567 524 // debug level 7
525 //
526
527 if(!fReconstructor->GetRecoParam() ){
528 AliError("Seed can not be used without a valid RecoParam.");
529 return kFALSE;
530 }
531
532 AliTRDchamberTimeBin *layer = 0x0;
a8276d32 533 if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7){
e8037fda 534 AliTRDtrackingChamber ch(*chamber);
535 ch.SetOwner();
29f95561 536 TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker);
537 cstreamer << "AttachClustersIter"
e8037fda 538 << "chamber.=" << &ch
29b87567 539 << "tracklet.=" << this
29b87567 540 << "\n";
541 }
542
35c24814 543 Float_t tquality;
29b87567 544 Double_t kroady = fReconstructor->GetRecoParam() ->GetRoad1y();
545 Double_t kroadz = fPadLength * .5 + 1.;
35c24814 546
547 // initialize configuration parameters
548 Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
549 Int_t niter = kZcorr ? 1 : 2;
550
29b87567 551 Double_t yexp, zexp;
552 Int_t ncl = 0;
35c24814 553 // start seed update
554 for (Int_t iter = 0; iter < niter; iter++) {
29b87567 555 ncl = 0;
556 for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
557 if(!(layer = chamber->GetTB(iTime))) continue;
558 if(!Int_t(*layer)) continue;
559
560 // define searching configuration
561 Double_t dxlayer = layer->GetX() - fX0;
562 if(c){
563 zexp = c->GetZ();
564 //Try 2 pad-rows in second iteration
565 if (iter > 0) {
566 zexp = fZref[0] + fZref[1] * dxlayer - zcorr;
567 if (zexp > c->GetZ()) zexp = c->GetZ() + fPadLength*0.5;
568 if (zexp < c->GetZ()) zexp = c->GetZ() - fPadLength*0.5;
569 }
570 } else zexp = fZref[0] + (kZcorr ? fZref[1] * dxlayer : 0.);
35c24814 571 yexp = fYref[0] + fYref[1] * dxlayer - zcorr;
29b87567 572
573 // Get and register cluster
574 Int_t index = layer->SearchNearestCluster(yexp, zexp, kroady, kroadz);
575 if (index < 0) continue;
576 AliTRDcluster *cl = (*layer)[index];
35c24814 577
29b87567 578 fIndexes[iTime] = layer->GetGlobalIndex(index);
579 fClusters[iTime] = cl;
580 fY[iTime] = cl->GetY();
581 fZ[iTime] = cl->GetZ();
582 ncl++;
583 }
35c24814 584 if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d ncl [%d] = %d", iter, fDet, ncl));
29b87567 585
586 if(ncl>1){
587 // calculate length of the time bin (calibration aware)
588 Int_t irp = 0; Float_t x[2]; Int_t tb[2];
589 for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
590 if(!fClusters[iTime]) continue;
591 x[irp] = fClusters[iTime]->GetX();
592 tb[irp] = iTime;
593 irp++;
594 if(irp==2) break;
595 }
596 fdX = (x[1] - x[0]) / (tb[0] - tb[1]);
597
598 // update X0 from the clusters (calibration/alignment aware)
599 for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
600 if(!(layer = chamber->GetTB(iTime))) continue;
601 if(!layer->IsT0()) continue;
602 if(fClusters[iTime]){
603 fX0 = fClusters[iTime]->GetX();
604 break;
605 } else { // we have to infere the position of the anode wire from the other clusters
606 for (Int_t jTime = iTime+1; jTime < AliTRDtrackerV1::GetNTimeBins(); jTime++) {
607 if(!fClusters[jTime]) continue;
608 fX0 = fClusters[jTime]->GetX() + fdX * (jTime - iTime);
f660dce9 609 break;
29b87567 610 }
29b87567 611 }
612 }
613
614 // update YZ reference point
615 // TODO
616
617 // update x reference positions (calibration/alignment aware)
618 for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
619 if(!fClusters[iTime]) continue;
0849f128 620 fX[iTime] = fX0 - fClusters[iTime]->GetX();
29b87567 621 }
622
623 AliTRDseed::Update();
624 }
35c24814 625 if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d nclFit [%d] = %d", iter, fDet, fN2));
29b87567 626
627 if(IsOK()){
628 tquality = GetQuality(kZcorr);
629 if(tquality < quality) break;
630 else quality = tquality;
631 }
632 kroadz *= 2.;
633 } // Loop: iter
634 if (!IsOK()) return kFALSE;
635
804bb02e 636 if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=1) CookLabels();
d937ad7a 637
b1957d3c 638 // set ExB angle
639 SetExB();
29b87567 640 UpdateUsed();
641 return kTRUE;
e4f2f73d 642}
643
644//____________________________________________________________________
b1957d3c 645Bool_t AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *chamber, Bool_t tilt)
e4f2f73d 646{
647 //
648 // Projective algorithm to attach clusters to seeding tracklets
649 //
650 // Parameters
651 //
652 // Output
653 //
654 // Detailed description
655 // 1. Collapse x coordinate for the full detector plane
656 // 2. truncated mean on y (r-phi) direction
657 // 3. purge clusters
658 // 4. truncated mean on z direction
659 // 5. purge clusters
660 // 6. fit tracklet
661 //
b1957d3c 662 Bool_t kPRINT = kFALSE;
29b87567 663 if(!fReconstructor->GetRecoParam() ){
664 AliError("Seed can not be used without a valid RecoParam.");
665 return kFALSE;
666 }
b1957d3c 667 // Initialize reco params for this tracklet
668 // 1. first time bin in the drift region
669 Int_t t0 = 4;
670 Int_t kClmin = Int_t(fReconstructor->GetRecoParam() ->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
29b87567 671
b1957d3c 672 Double_t syRef = TMath::Sqrt(fRefCov[0]);
29b87567 673 //define roads
b1957d3c 674 Double_t kroady = 1.;
675 //fReconstructor->GetRecoParam() ->GetRoad1y();
29b87567 676 Double_t kroadz = fPadLength * 1.5 + 1.;
b1957d3c 677 if(kPRINT) printf("AttachClusters() sy[%f] road[%f]\n", syRef, kroady);
29b87567 678
679 // working variables
b1957d3c 680 const Int_t kNrows = 16;
681 AliTRDcluster *clst[kNrows][knTimebins];
682 Double_t cond[4], dx, dy, yt, zt,
683 yres[kNrows][knTimebins];
684 Int_t idxs[kNrows][knTimebins], ncl[kNrows], ncls = 0;
685 memset(ncl, 0, kNrows*sizeof(Int_t));
686 memset(clst, 0, kNrows*knTimebins*sizeof(AliTRDcluster*));
687
29b87567 688 // Do cluster projection
b1957d3c 689 AliTRDcluster *c = 0x0;
29b87567 690 AliTRDchamberTimeBin *layer = 0x0;
b1957d3c 691 Bool_t kBUFFER = kFALSE;
692 for (Int_t it = 0; it < AliTRDtrackerV1::GetNTimeBins(); it++) {
693 if(!(layer = chamber->GetTB(it))) continue;
29b87567 694 if(!Int_t(*layer)) continue;
695
b1957d3c 696 dx = fX0 - layer->GetX();
697 yt = fYref[0] - fYref[1] * dx;
698 zt = fZref[0] - fZref[1] * dx;
699 if(kPRINT) printf("\t%2d dx[%f] yt[%f] zt[%f]\n", it, dx, yt, zt);
700
701 // select clusters on a 5 sigmaKalman level
702 cond[0] = yt; cond[2] = kroady;
703 cond[1] = zt; cond[3] = kroadz;
704 Int_t n=0, idx[6];
705 layer->GetClusters(cond, idx, n, 6);
706 for(Int_t ic = n; ic--;){
707 c = (*layer)[idx[ic]];
708 dy = yt - c->GetY();
709 dy += tilt ? fTilt * (c->GetZ() - zt) : 0.;
710 // select clusters on a 3 sigmaKalman level
711/* if(tilt && TMath::Abs(dy) > 3.*syRef){
712 printf("too large !!!\n");
713 continue;
714 }*/
715 Int_t r = c->GetPadRow();
716 if(kPRINT) printf("\t\t%d dy[%f] yc[%f] r[%d]\n", ic, TMath::Abs(dy), c->GetY(), r);
717 clst[r][ncl[r]] = c;
718 idxs[r][ncl[r]] = idx[ic];
719 yres[r][ncl[r]] = dy;
720 ncl[r]++; ncls++;
721
722 if(ncl[r] >= knTimebins) {
723 AliWarning(Form("Cluster candidates reached limit %d. Some may be lost.", knTimebins));
724 kBUFFER = kTRUE;
29b87567 725 break;
726 }
727 }
b1957d3c 728 if(kBUFFER) break;
29b87567 729 }
b1957d3c 730 if(kPRINT) printf("Found %d clusters\n", ncls);
731 if(ncls<kClmin) return kFALSE;
732
733 // analyze each row individualy
734 Double_t mean, syDis;
735 Int_t nrow[] = {0, 0, 0}, nr = 0, lr=-1;
736 for(Int_t ir=kNrows; ir--;){
737 if(!(ncl[ir])) continue;
738 if(lr>0 && lr-ir != 1){
739 if(kPRINT) printf("W - gap in rows attached !!\n");
29b87567 740 }
b1957d3c 741 if(kPRINT) printf("\tir[%d] lr[%d] n[%d]\n", ir, lr, ncl[ir]);
742 // Evaluate truncated mean on the y direction
743 if(ncl[ir] > 3) AliMathBase::EvaluateUni(ncl[ir], yres[ir], mean, syDis, Int_t(ncl[ir]*.8));
744 else {
745 mean = 0.; syDis = 0.;
746 }
747
748 // TODO check mean and sigma agains cluster resolution !!
749 if(kPRINT) printf("\tr[%2d] m[%f %5.3fsigma] s[%f]\n", ir, mean, TMath::Abs(mean/syRef), syDis);
750 // select clusters on a 3 sigmaDistr level
751 Bool_t kFOUND = kFALSE;
752 for(Int_t ic = ncl[ir]; ic--;){
753 if(yres[ir][ic] - mean > 3. * syDis){
754 clst[ir][ic] = 0x0; continue;
755 }
756 nrow[nr]++; kFOUND = kTRUE;
757 }
758 // exit loop
759 if(kFOUND) nr++;
760 lr = ir; if(nr>=3) break;
29b87567 761 }
b1957d3c 762 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]);
763
764 // classify cluster rows
765 Int_t row = -1;
766 switch(nr){
767 case 1:
768 row = lr;
769 break;
770 case 2:
771 SetBit(kRowCross, kTRUE); // mark pad row crossing
772 if(nrow[0] > nrow[1]){ row = lr+1; lr = -1;}
773 else{
774 row = lr; lr = 1;
775 nrow[2] = nrow[1];
776 nrow[1] = nrow[0];
777 nrow[0] = nrow[2];
29b87567 778 }
b1957d3c 779 break;
780 case 3:
781 SetBit(kRowCross, kTRUE); // mark pad row crossing
782 break;
29b87567 783 }
b1957d3c 784 if(kPRINT) printf("\trow[%d] n[%d]\n\n", row, nrow[0]);
785 if(row<0) return kFALSE;
29b87567 786
b1957d3c 787 // Select and store clusters
788 // We should consider here :
789 // 1. How far is the chamber boundary
790 // 2. How big is the mean
29b87567 791 fN2 = 0;
b1957d3c 792 for (Int_t ir = 0; ir < nr; ir++) {
793 Int_t jr = row + ir*lr;
794 if(kPRINT) printf("\tattach %d clusters for row %d\n", ncl[jr], jr);
795 for (Int_t ic = 0; ic < ncl[jr]; ic++) {
796 if(!(c = clst[jr][ic])) continue;
797 Int_t it = c->GetPadTime();
798 // TODO proper indexing of clusters !!
799 fIndexes[it+35*ir] = chamber->GetTB(it)->GetGlobalIndex(idxs[jr][ic]);
800 fClusters[it+35*ir] = c;
29b87567 801
b1957d3c 802 //printf("\tid[%2d] it[%d] idx[%d]\n", ic, it, fIndexes[it]);
803
804 fN2++;
805 }
806 }
807
29b87567 808 // number of minimum numbers of clusters expected for the tracklet
e4f2f73d 809 if (fN2 < kClmin){
29b87567 810 AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin));
e4f2f73d 811 fN2 = 0;
812 return kFALSE;
813 }
0906e73e 814
b1957d3c 815 // update used clusters and select
29b87567 816 fNUsed = 0;
b1957d3c 817 for (Int_t it = 0; it < AliTRDtrackerV1::GetNTimeBins(); it++) {
818 if(fClusters[it] && fClusters[it]->IsUsed()) fNUsed++;
819 if(fClusters[it+35] && fClusters[it+35]->IsUsed()) fNUsed++;
29b87567 820 }
0906e73e 821 if (fN2-fNUsed < kClmin){
b1957d3c 822 //AliWarning(Form("Too many clusters already in use %d (from %d).", fNUsed, fN2));
0906e73e 823 fN2 = 0;
824 return kFALSE;
825 }
b1957d3c 826
827 // set the Lorentz angle for this tracklet
828 SetExB();
829
830 // calculate dx for time bins in the drift region (calibration aware)
d86ed84c 831 Int_t irp = 0; Float_t x[2]={0., 0.}; Int_t tb[2] = {0, 0};
b1957d3c 832 for (Int_t it = t0; it < AliTRDtrackerV1::GetNTimeBins(); it++) {
833 if(!fClusters[it]) continue;
834 x[irp] = fClusters[it]->GetX();
835 tb[irp] = it;
836 irp++;
837 if(irp==2) break;
d86ed84c 838 }
839 Int_t dtb = tb[1] - tb[0];
840 fdX = dtb ? (x[0] - x[1]) / dtb : 0.15;
b1957d3c 841
842 // update X0 from the clusters (calibration/alignment aware) TODO remove dependence on x0 !!
843 for (Int_t it = 0; it < AliTRDtrackerV1::GetNTimeBins(); it++) {
844 if(!(layer = chamber->GetTB(it))) continue;
845 if(!layer->IsT0()) continue;
846 if(fClusters[it]){
847 fX0 = fClusters[it]->GetX();
848 break;
849 } else { // we have to infere the position of the anode wire from the other clusters
850 for (Int_t jt = it+1; jt < AliTRDtrackerV1::GetNTimeBins(); jt++) {
851 if(!fClusters[jt]) continue;
852 fX0 = fClusters[jt]->GetX() + fdX * (jt - it);
853 break;
854 }
855 }
856 }
857
29b87567 858 return kTRUE;
e4f2f73d 859}
860
03cef9b2 861//____________________________________________________________
862void AliTRDseedV1::Bootstrap(const AliTRDReconstructor *rec)
863{
864// Fill in all derived information. It has to be called after recovery from file or HLT.
865// The primitive data are
866// - list of clusters
867// - detector (as the detector will be removed from clusters)
868// - position of anode wire (fX0) - temporary
869// - track reference position and direction
870// - momentum of the track
871// - time bin length [cm]
872//
873// A.Bercuci <A.Bercuci@gsi.de> Oct 30th 2008
874//
875 fReconstructor = rec;
876 AliTRDgeometry g;
877 AliTRDpadPlane *pp = g.GetPadPlane(fDet);
878 fTilt = TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle());
879 fPadLength = pp->GetLengthIPad();
880 fSnp = fYref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);
881 fTgl = fZref[1];
882 fN = 0; fN2 = 0; fMPads = 0.;
883 AliTRDcluster **cit = &fClusters[0];
884 for(Int_t ic = knTimebins; ic--; cit++){
885 if(!(*cit)) return;
886 fN++; fN2++;
887 fX[ic] = (*cit)->GetX() - fX0;
888 fY[ic] = (*cit)->GetY();
889 fZ[ic] = (*cit)->GetZ();
890 }
891 Update(); // Fit();
892 CookLabels();
893 GetProbability();
894}
895
896
e4f2f73d 897//____________________________________________________________________
d937ad7a 898Bool_t AliTRDseedV1::Fit(Bool_t tilt, Int_t errors)
e4f2f73d 899{
900 //
901 // Linear fit of the tracklet
902 //
903 // Parameters :
904 //
905 // Output :
906 // True if successful
907 //
908 // Detailed description
909 // 2. Check if tracklet crosses pad row boundary
910 // 1. Calculate residuals in the y (r-phi) direction
911 // 3. Do a Least Square Fit to the data
912 //
913
29b87567 914 const Int_t kClmin = 8;
010d62b0 915
9462866a 916
917 // cluster error parametrization parameters
010d62b0 918 // 1. sy total charge
9462866a 919 const Float_t sq0inv = 0.019962; // [1/q0]
920 const Float_t sqb = 1.0281564; //[cm]
010d62b0 921 // 2. sy for the PRF
922 const Float_t scy[AliTRDgeometry::kNlayer][4] = {
d937ad7a 923 {2.827e-02, 9.600e-04, 4.296e-01, 2.271e-02},
924 {2.952e-02,-2.198e-04, 4.146e-01, 2.339e-02},
925 {3.090e-02, 1.514e-03, 4.020e-01, 2.402e-02},
926 {3.260e-02,-2.037e-03, 3.946e-01, 2.509e-02},
927 {3.439e-02,-3.601e-04, 3.883e-01, 2.623e-02},
928 {3.510e-02, 2.066e-03, 3.651e-01, 2.588e-02},
010d62b0 929 };
930 // 3. sy parallel to the track
d937ad7a 931 const Float_t sy0 = 2.649e-02; // [cm]
932 const Float_t sya = -8.864e-04; // [cm]
933 const Float_t syb = -2.435e-01; // [cm]
934
010d62b0 935 // 4. sx parallel to the track
d937ad7a 936 const Float_t sxgc = 5.427e-02;
937 const Float_t sxgm = 7.783e-01;
938 const Float_t sxgs = 2.743e-01;
939 const Float_t sxe0 =-2.065e+00;
940 const Float_t sxe1 =-2.978e-02;
941
010d62b0 942 // 5. sx perpendicular to the track
d937ad7a 943// const Float_t sxd0 = 1.881e-02;
944// const Float_t sxd1 =-4.101e-01;
945// const Float_t sxd2 = 1.572e+00;
946
2f7d6ac8 947 // get track direction
948 Double_t y0 = fYref[0];
949 Double_t dydx = fYref[1];
950 Double_t z0 = fZref[0];
951 Double_t dzdx = fZref[1];
952 Double_t yt, zt;
ae4e8b84 953
29b87567 954 const Int_t kNtb = AliTRDtrackerV1::GetNTimeBins();
b1957d3c 955 //AliTRDtrackerV1::AliTRDLeastSquare fitterZ;
24d8660e 956 TLinearFitter fitterY(1, "pol1");
29b87567 957 // convertion factor from square to gauss distribution for sigma
b1957d3c 958 //Double_t convert = 1./TMath::Sqrt(12.);
ae4e8b84 959
29b87567 960 // book cluster information
b1957d3c 961 Double_t q, xc[knTimebins], yc[knTimebins], zc[knTimebins], sy[knTimebins]/*, sz[knTimebins]*/;
962// Int_t zRow[knTimebins];
9462866a 963
010d62b0 964 Int_t ily = AliTRDgeometry::GetLayer(fDet);
b1957d3c 965 fN = 0; //fXref = 0.; Double_t ssx = 0.;
9eb2d46c 966 AliTRDcluster *c=0x0, **jc = &fClusters[0];
9eb2d46c 967 for (Int_t ic=0; ic<kNtb; ic++, ++jc) {
b1957d3c 968 //zRow[ic] = -1;
29b87567 969 xc[ic] = -1.;
970 yc[ic] = 999.;
971 zc[ic] = 999.;
972 sy[ic] = 0.;
b1957d3c 973 //sz[ic] = 0.;
9eb2d46c 974 if(!(c = (*jc))) continue;
29b87567 975 if(!c->IsInChamber()) continue;
9462866a 976
29b87567 977 Float_t w = 1.;
978 if(c->GetNPads()>4) w = .5;
979 if(c->GetNPads()>5) w = .2;
010d62b0 980
b1957d3c 981 //zRow[fN] = c->GetPadRow();
d937ad7a 982 // correct cluster position for PRF and v drift
983 Double_t x, y; GetClusterXY(c, x, y);
984 xc[fN] = fX0 - x;
985 yc[fN] = y;
2f7d6ac8 986 zc[fN] = c->GetZ();
987
988 // extrapolated y value for the track
989 yt = y0 - xc[fN]*dydx;
990 // extrapolated z value for the track
991 zt = z0 - xc[fN]*dzdx;
992 // tilt correction
993 if(tilt) yc[fN] -= fTilt*(zc[fN] - zt);
994
010d62b0 995 // ELABORATE CLUSTER ERROR
996 // TODO to be moved to AliTRDcluster
9462866a 997 q = TMath::Abs(c->GetQ());
d937ad7a 998 Double_t tgg = (dydx-fExB)/(1.+dydx*fExB); tgg *= tgg;
010d62b0 999 // basic y error (|| to track).
d937ad7a 1000 sy[fN] = xc[fN] < AliTRDgeometry::CamHght() ? 2. : sy0 + sya*TMath::Exp(1./(xc[fN]+syb));
1001 //printf("cluster[%d]\n\tsy[0] = %5.3e [um]\n", fN, sy[fN]*1.e4);
010d62b0 1002 // y error due to total charge
1003 sy[fN] += sqb*(1./q - sq0inv);
d937ad7a 1004 //printf("\tsy[1] = %5.3e [um]\n", sy[fN]*1.e4);
010d62b0 1005 // y error due to PRF
1006 sy[fN] += scy[ily][0]*TMath::Gaus(c->GetCenter(), scy[ily][1], scy[ily][2]) - scy[ily][3];
d937ad7a 1007 //printf("\tsy[2] = %5.3e [um]\n", sy[fN]*1.e4);
1008
010d62b0 1009 sy[fN] *= sy[fN];
1010
1011 // ADD ERROR ON x
9462866a 1012 // error of drift length parallel to the track
d937ad7a 1013 Double_t sx = sxgc*TMath::Gaus(xc[fN], sxgm, sxgs) + TMath::Exp(sxe0+sxe1*xc[fN]); // [cm]
1014 //printf("\tsx[0] = %5.3e [um]\n", sx*1.e4);
9462866a 1015 // error of drift length perpendicular to the track
1016 //sx += sxd0 + sxd1*d + sxd2*d*d;
d937ad7a 1017 sx *= sx; // square sx
1018 // update xref
b1957d3c 1019 //fXref += xc[fN]/sx; ssx+=1./sx;
d937ad7a 1020
9462866a 1021 // add error from ExB
d937ad7a 1022 if(errors>0) sy[fN] += fExB*fExB*sx;
1023 //printf("\tsy[3] = %5.3e [um^2]\n", sy[fN]*1.e8);
1024
1025 // global radial error due to misalignment/miscalibration
1026 Double_t sx0 = 0.; sx0 *= sx0;
1027 // add sx contribution to sy due to track angle
1028 if(errors>1) sy[fN] += tgg*(sx+sx0);
1029 // TODO we should add tilt pad correction here
1030 //printf("\tsy[4] = %5.3e [um^2]\n", sy[fN]*1.e8);
1031 c->SetSigmaY2(sy[fN]);
1032
9462866a 1033 sy[fN] = TMath::Sqrt(sy[fN]);
24d8660e 1034 fitterY.AddPoint(&xc[fN], yc[fN]/*-yt*/, sy[fN]);
2f7d6ac8 1035 fN++;
29b87567 1036 }
47d5d320 1037 // to few clusters
2f7d6ac8 1038 if (fN < kClmin) return kFALSE;
1039
d937ad7a 1040 // fit XY
2f7d6ac8 1041 fitterY.Eval();
d937ad7a 1042 fYfit[0] = fitterY.GetParameter(0);
1043 fYfit[1] = -fitterY.GetParameter(1);
1044 // store covariance
1045 Double_t *p = fitterY.GetCovarianceMatrix();
1046 fCov[0] = p[0]; // variance of y0
1047 fCov[1] = p[1]; // covariance of y0, dydx
1048 fCov[2] = p[3]; // variance of dydx
b1957d3c 1049 // the ref radial position is set at the minimum of
1050 // the y variance of the tracklet
1051 fXref = -fCov[1]/fCov[2]; //fXref = fX0 - fXref;
1052
1053 // fit XZ
1054 if(IsRowCross()){
1055 // TODO pad row cross position estimation !!!
1056 //AliInfo(Form("Padrow cross in detector %d", fDet));
1057 fZfit[0] = .5*(zc[0]+zc[fN-1]); fZfit[1] = 0.;
1058 } else {
1059 fZfit[0] = zc[0]; fZfit[1] = 0.;
29b87567 1060 }
1061
29b87567 1062
b1957d3c 1063// // determine z offset of the fit
1064// Float_t zslope = 0.;
1065// Int_t nchanges = 0, nCross = 0;
1066// if(nz==2){ // tracklet is crossing pad row
1067// // Find the break time allowing one chage on pad-rows
1068// // with maximal number of accepted clusters
1069// Int_t padRef = zRow[0];
1070// for (Int_t ic=1; ic<fN; ic++) {
1071// if(zRow[ic] == padRef) continue;
1072//
1073// // debug
1074// if(zRow[ic-1] == zRow[ic]){
1075// printf("ERROR in pad row change!!!\n");
1076// }
1077//
1078// // evaluate parameters of the crossing point
1079// Float_t sx = (xc[ic-1] - xc[ic])*convert;
1080// fCross[0] = .5 * (xc[ic-1] + xc[ic]);
1081// fCross[2] = .5 * (zc[ic-1] + zc[ic]);
1082// fCross[3] = TMath::Max(dzdx * sx, .01);
1083// zslope = zc[ic-1] > zc[ic] ? 1. : -1.;
1084// padRef = zRow[ic];
1085// nCross = ic;
1086// nchanges++;
1087// }
1088// }
1089//
1090// // condition on nCross and reset nchanges TODO
1091//
1092// if(nchanges==1){
1093// if(dzdx * zslope < 0.){
1094// AliInfo("Tracklet-Track mismatch in dzdx. TODO.");
1095// }
1096//
1097//
1098// //zc[nc] = fitterZ.GetFunctionParameter(0);
1099// fCross[1] = fYfit[0] - fCross[0] * fYfit[1];
1100// fCross[0] = fX0 - fCross[0];
1101// }
29b87567 1102
2389e96f 1103 UpdateUsed();
29b87567 1104 return kTRUE;
e4f2f73d 1105}
1106
e4f2f73d 1107
1108//___________________________________________________________________
203967fc 1109void AliTRDseedV1::Print(Option_t *o) const
e4f2f73d 1110{
1111 //
1112 // Printing the seedstatus
1113 //
1114
203967fc 1115 AliInfo(Form("Det[%3d] Tilt[%+6.2f] Pad[%5.2f]", fDet, fTilt, fPadLength));
1116 AliInfo(Form("Nattach[%2d] Nfit[%2d] Nuse[%2d] pads[%f]", fN, fN2, fNUsed, fMPads));
1117 AliInfo(Form("x[%7.2f] y[%7.2f] z[%7.2f] dydx[%5.2f] dzdx[%5.2f]", fX0, fYfit[0], fZfit[0], fYfit[1], fZfit[1]));
1118 AliInfo(Form("Ref y[%7.2f] z[%7.2f] dydx[%5.2f] dzdx[%5.2f]", fYref[0], fZref[0], fYref[1], fZref[1]))
1119
1120
1121 if(strcmp(o, "a")!=0) return;
1122
4dc4dc2e 1123 AliTRDcluster* const* jc = &fClusters[0];
1124 for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++, jc++) {
1125 if(!(*jc)) continue;
203967fc 1126 (*jc)->Print(o);
4dc4dc2e 1127 }
e4f2f73d 1128}
47d5d320 1129
203967fc 1130
1131//___________________________________________________________________
1132Bool_t AliTRDseedV1::IsEqual(const TObject *o) const
1133{
1134 // Checks if current instance of the class has the same essential members
1135 // as the given one
1136
1137 if(!o) return kFALSE;
1138 const AliTRDseedV1 *inTracklet = dynamic_cast<const AliTRDseedV1*>(o);
1139 if(!inTracklet) return kFALSE;
1140
1141 for (Int_t i = 0; i < 2; i++){
1142 if ( fYref[i] != inTracklet->GetYref(i) ) return kFALSE;
1143 if ( fZref[i] != inTracklet->GetZref(i) ) return kFALSE;
1144 }
1145
1146 if ( fSigmaY != inTracklet->GetSigmaY() ) return kFALSE;
1147 if ( fSigmaY2 != inTracklet->GetSigmaY2() ) return kFALSE;
1148 if ( fTilt != inTracklet->GetTilt() ) return kFALSE;
1149 if ( fPadLength != inTracklet->GetPadLength() ) return kFALSE;
1150
1151 for (Int_t i = 0; i < knTimebins; i++){
1152 if ( fX[i] != inTracklet->GetX(i) ) return kFALSE;
1153 if ( fY[i] != inTracklet->GetY(i) ) return kFALSE;
1154 if ( fZ[i] != inTracklet->GetZ(i) ) return kFALSE;
1155 if ( fIndexes[i] != inTracklet->GetIndexes(i) ) return kFALSE;
1156 if ( fUsable[i] != inTracklet->IsUsable(i) ) return kFALSE;
1157 }
1158
1159 for (Int_t i=0; i < 2; i++){
1160 if ( fYfit[i] != inTracklet->GetYfit(i) ) return kFALSE;
1161 if ( fZfit[i] != inTracklet->GetZfit(i) ) return kFALSE;
1162 if ( fYfitR[i] != inTracklet->GetYfitR(i) ) return kFALSE;
1163 if ( fZfitR[i] != inTracklet->GetZfitR(i) ) return kFALSE;
1164 if ( fLabels[i] != inTracklet->GetLabels(i) ) return kFALSE;
1165 }
1166
1167 if ( fMeanz != inTracklet->GetMeanz() ) return kFALSE;
1168 if ( fZProb != inTracklet->GetZProb() ) return kFALSE;
1169 if ( fN2 != inTracklet->GetN2() ) return kFALSE;
1170 if ( fNUsed != inTracklet->GetNUsed() ) return kFALSE;
1171 if ( fFreq != inTracklet->GetFreq() ) return kFALSE;
1172 if ( fNChange != inTracklet->GetNChange() ) return kFALSE;
1173 if ( fNChange != inTracklet->GetNChange() ) return kFALSE;
1174
1175 if ( fC != inTracklet->GetC() ) return kFALSE;
1176 if ( fCC != inTracklet->GetCC() ) return kFALSE;
1177 if ( fChi2 != inTracklet->GetChi2() ) return kFALSE;
1178 // if ( fChi2Z != inTracklet->GetChi2Z() ) return kFALSE;
1179
1180 if ( fDet != inTracklet->GetDetector() ) return kFALSE;
1181 if ( fMom != inTracklet->GetMomentum() ) return kFALSE;
1182 if ( fdX != inTracklet->GetdX() ) return kFALSE;
1183
1184 for (Int_t iCluster = 0; iCluster < knTimebins; iCluster++){
1185 AliTRDcluster *curCluster = fClusters[iCluster];
1186 AliTRDcluster *inCluster = inTracklet->GetClusters(iCluster);
1187 if (curCluster && inCluster){
1188 if (! curCluster->IsEqual(inCluster) ) {
1189 curCluster->Print();
1190 inCluster->Print();
1191 return kFALSE;
1192 }
1193 } else {
1194 // if one cluster exists, and corresponding
1195 // in other tracklet doesn't - return kFALSE
1196 if(curCluster || inCluster) return kFALSE;
1197 }
1198 }
1199 return kTRUE;
1200}