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