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