dding default trigger description and trigger domain parameters for AliHLTTriggerMenu...
[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)
e3cf3d02 58 :TObject()
3a039a31 59 ,fReconstructor(0x0)
ae4e8b84 60 ,fClusterIter(0x0)
e3cf3d02 61 ,fExB(0.)
62 ,fVD(0.)
63 ,fT0(0.)
64 ,fS2PRF(0.)
65 ,fDiffL(0.)
66 ,fDiffT(0.)
ae4e8b84 67 ,fClusterIdx(0)
f29f13a6 68// ,fUsable(0)
e3cf3d02 69 ,fN2(0)
70 ,fNUsed(0)
ae4e8b84 71 ,fDet(det)
e3cf3d02 72 ,fTilt(0.)
73 ,fPadLength(0.)
0906e73e 74 ,fMom(0.)
bcb6fb78 75 ,fdX(0.)
e3cf3d02 76 ,fX0(0.)
77 ,fX(0.)
78 ,fY(0.)
79 ,fZ(0.)
80 ,fS2Y(0.)
81 ,fS2Z(0.)
82 ,fC(0.)
83 ,fChi2(0.)
e4f2f73d 84{
85 //
86 // Constructor
87 //
e3cf3d02 88 for(Int_t ic=kNTimeBins; ic--;) fIndexes[ic] = -1;
89 memset(fClusters, 0, kNTimeBins*sizeof(AliTRDcluster*));
90 fYref[0] = 0.; fYref[1] = 0.;
91 fZref[0] = 0.; fZref[1] = 0.;
92 fYfit[0] = 0.; fYfit[1] = 0.;
93 fZfit[0] = 0.; fZfit[1] = 0.;
94 memset(fdEdx, 0, kNSlices*sizeof(Float_t));
29b87567 95 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = -1.;
e3cf3d02 96 fLabels[0]=-1; fLabels[1]=-1; // most freq MC labels
97 fLabels[2]=0; // number of different labels for tracklet
6e4d4425 98 fRefCov[0] = 1.; fRefCov[1] = 0.; fRefCov[2] = 1.;
d937ad7a 99 // covariance matrix [diagonal]
100 // default sy = 200um and sz = 2.3 cm
101 fCov[0] = 4.e-4; fCov[1] = 0.; fCov[2] = 5.3;
f29f13a6 102 SetStandAlone(kFALSE);
e4f2f73d 103}
104
105//____________________________________________________________________
0906e73e 106AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
e3cf3d02 107 :TObject((TObject&)ref)
108 ,fReconstructor(0x0)
ae4e8b84 109 ,fClusterIter(0x0)
e3cf3d02 110 ,fExB(0.)
111 ,fVD(0.)
112 ,fT0(0.)
113 ,fS2PRF(0.)
114 ,fDiffL(0.)
115 ,fDiffT(0.)
ae4e8b84 116 ,fClusterIdx(0)
f29f13a6 117// ,fUsable(0)
e3cf3d02 118 ,fN2(0)
119 ,fNUsed(0)
120 ,fDet(-1)
121 ,fTilt(0.)
122 ,fPadLength(0.)
123 ,fMom(0.)
124 ,fdX(0.)
125 ,fX0(0.)
126 ,fX(0.)
127 ,fY(0.)
128 ,fZ(0.)
129 ,fS2Y(0.)
130 ,fS2Z(0.)
131 ,fC(0.)
132 ,fChi2(0.)
e4f2f73d 133{
134 //
135 // Copy Constructor performing a deep copy
136 //
e3cf3d02 137 if(this != &ref){
138 ref.Copy(*this);
139 }
29b87567 140 SetBit(kOwner, kFALSE);
f29f13a6 141 SetStandAlone(ref.IsStandAlone());
fbb2ea06 142}
d9950a5a 143
0906e73e 144
e4f2f73d 145//____________________________________________________________________
146AliTRDseedV1& AliTRDseedV1::operator=(const AliTRDseedV1 &ref)
147{
148 //
149 // Assignment Operator using the copy function
150 //
151
29b87567 152 if(this != &ref){
153 ref.Copy(*this);
154 }
221ab7e0 155 SetBit(kOwner, kFALSE);
156
29b87567 157 return *this;
e4f2f73d 158}
159
160//____________________________________________________________________
161AliTRDseedV1::~AliTRDseedV1()
162{
163 //
164 // Destructor. The RecoParam object belongs to the underlying tracker.
165 //
166
29b87567 167 //printf("I-AliTRDseedV1::~AliTRDseedV1() : Owner[%s]\n", IsOwner()?"YES":"NO");
e4f2f73d 168
e3cf3d02 169 if(IsOwner()) {
170 for(int itb=0; itb<kNTimeBins; itb++){
29b87567 171 if(!fClusters[itb]) continue;
172 //AliInfo(Form("deleting c %p @ %d", fClusters[itb], itb));
173 delete fClusters[itb];
174 fClusters[itb] = 0x0;
175 }
e3cf3d02 176 }
e4f2f73d 177}
178
179//____________________________________________________________________
180void AliTRDseedV1::Copy(TObject &ref) const
181{
182 //
183 // Copy function
184 //
185
29b87567 186 //AliInfo("");
187 AliTRDseedV1 &target = (AliTRDseedV1 &)ref;
188
e3cf3d02 189 target.fReconstructor = fReconstructor;
ae4e8b84 190 target.fClusterIter = 0x0;
e3cf3d02 191 target.fExB = fExB;
192 target.fVD = fVD;
193 target.fT0 = fT0;
194 target.fS2PRF = fS2PRF;
195 target.fDiffL = fDiffL;
196 target.fDiffT = fDiffT;
ae4e8b84 197 target.fClusterIdx = 0;
f29f13a6 198// target.fUsable = fUsable;
e3cf3d02 199 target.fN2 = fN2;
200 target.fNUsed = fNUsed;
ae4e8b84 201 target.fDet = fDet;
e3cf3d02 202 target.fTilt = fTilt;
203 target.fPadLength = fPadLength;
29b87567 204 target.fMom = fMom;
29b87567 205 target.fdX = fdX;
e3cf3d02 206 target.fX0 = fX0;
207 target.fX = fX;
208 target.fY = fY;
209 target.fZ = fZ;
210 target.fS2Y = fS2Y;
211 target.fS2Z = fS2Z;
212 target.fC = fC;
213 target.fChi2 = fChi2;
29b87567 214
e3cf3d02 215 memcpy(target.fIndexes, fIndexes, kNTimeBins*sizeof(Int_t));
216 memcpy(target.fClusters, fClusters, kNTimeBins*sizeof(AliTRDcluster*));
217 target.fYref[0] = fYref[0]; target.fYref[1] = fYref[1];
218 target.fZref[0] = fZref[0]; target.fZref[1] = fZref[1];
219 target.fYfit[0] = fYfit[0]; target.fYfit[1] = fYfit[1];
220 target.fZfit[0] = fZfit[0]; target.fZfit[1] = fZfit[1];
221 memcpy(target.fdEdx, fdEdx, kNSlices*sizeof(Float_t));
222 memcpy(target.fProb, fProb, AliPID::kSPECIES*sizeof(Float_t));
223 memcpy(target.fLabels, fLabels, 3*sizeof(Int_t));
224 memcpy(target.fRefCov, fRefCov, 3*sizeof(Double_t));
225 memcpy(target.fCov, fCov, 3*sizeof(Double_t));
29b87567 226
e3cf3d02 227 TObject::Copy(ref);
e4f2f73d 228}
229
0906e73e 230
231//____________________________________________________________
f3d3af1b 232Bool_t AliTRDseedV1::Init(AliTRDtrackV1 *track)
0906e73e 233{
234// Initialize this tracklet using the track information
235//
236// Parameters:
237// track - the TRD track used to initialize the tracklet
238//
239// Detailed description
240// The function sets the starting point and direction of the
241// tracklet according to the information from the TRD track.
242//
243// Caution
244// The TRD track has to be propagated to the beginning of the
245// chamber where the tracklet will be constructed
246//
247
29b87567 248 Double_t y, z;
249 if(!track->GetProlongation(fX0, y, z)) return kFALSE;
b1957d3c 250 UpDate(track);
29b87567 251 return kTRUE;
0906e73e 252}
253
bcb6fb78 254
e3cf3d02 255//_____________________________________________________________________________
256void AliTRDseedV1::Reset()
257{
258 //
259 // Reset seed
260 //
261 fExB=0.;fVD=0.;fT0=0.;fS2PRF=0.;
262 fDiffL=0.;fDiffT=0.;
f29f13a6 263 fClusterIdx=0;//fUsable=0;
264 fN2=0; fNUsed=0;
e3cf3d02 265 fDet=-1;fTilt=0.;fPadLength=0.;
266 fMom=0.;
267 fdX=0.;fX0=0.; fX=0.; fY=0.; fZ=0.;
268 fS2Y=0.; fS2Z=0.;
269 fC=0.; fChi2 = 0.;
270
271 for(Int_t ic=kNTimeBins; ic--;) fIndexes[ic] = -1;
272 memset(fClusters, 0, kNTimeBins*sizeof(AliTRDcluster*));
273 fYref[0] = 0.; fYref[1] = 0.;
274 fZref[0] = 0.; fZref[1] = 0.;
275 fYfit[0] = 0.; fYfit[1] = 0.;
276 fZfit[0] = 0.; fZfit[1] = 0.;
277 memset(fdEdx, 0, kNSlices*sizeof(Float_t));
278 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = -1.;
279 fLabels[0]=-1; fLabels[1]=-1; // most freq MC labels
280 fLabels[2]=0; // number of different labels for tracklet
281 fRefCov[0] = 1.; fRefCov[1] = 0.; fRefCov[2] = 1.;
282 // covariance matrix [diagonal]
283 // default sy = 200um and sz = 2.3 cm
284 fCov[0] = 4.e-4; fCov[1] = 0.; fCov[2] = 5.3;
285}
286
bcb6fb78 287//____________________________________________________________________
b1957d3c 288void AliTRDseedV1::UpDate(const AliTRDtrackV1 *trk)
289{
290 // update tracklet reference position from the TRD track
291 // Funny name to avoid the clash with the function AliTRDseed::Update() (has to be made obselete)
292
e3cf3d02 293 Double_t fSnp = trk->GetSnp();
294 Double_t fTgl = trk->GetTgl();
b1957d3c 295 fMom = trk->GetP();
296 fYref[1] = fSnp/(1. - fSnp*fSnp);
297 fZref[1] = fTgl;
298 SetCovRef(trk->GetCovariance());
299
300 Double_t dx = trk->GetX() - fX0;
301 fYref[0] = trk->GetY() - dx*fYref[1];
302 fZref[0] = trk->GetZ() - dx*fZref[1];
303}
304
e3cf3d02 305//_____________________________________________________________________________
306void AliTRDseedV1::UpdateUsed()
307{
308 //
f29f13a6 309 // Calculate number of used clusers in the tracklet
e3cf3d02 310 //
311
312 fNUsed = 0;
313 for (Int_t i = kNTimeBins; i--; ) {
314 if (!fClusters[i]) continue;
f29f13a6 315 if(fClusters[i]->IsUsed()) fNUsed++;
316 else if(fClusters[i]->IsShared() && IsStandAlone()) fNUsed++;
e3cf3d02 317 }
318}
319
320//_____________________________________________________________________________
321void AliTRDseedV1::UseClusters()
322{
323 //
324 // Use clusters
325 //
f29f13a6 326 // In stand alone mode:
327 // Clusters which are marked as used or shared from another track are
328 // removed from the tracklet
329 //
330 // In barrel mode:
331 // - Clusters which are used by another track become shared
332 // - Clusters which are attached to a kink track become shared
333 //
e3cf3d02 334 AliTRDcluster **c = &fClusters[0];
335 for (Int_t ic=kNTimeBins; ic--; c++) {
336 if(!(*c)) continue;
f29f13a6 337 if(IsStandAlone()){
338 if((*c)->IsShared() || (*c)->IsUsed()){
339 (*c) = 0x0;
340 fIndexes[ic] = -1;
341 fN2--;
342 continue;
343 }
344 }
345 else{
346 if((*c)->IsUsed() || IsKink()){
347 (*c)->SetShared();
348 continue;
349 }
350 }
351 (*c)->Use();
e3cf3d02 352 }
353}
354
355
f29f13a6 356
b1957d3c 357//____________________________________________________________________
bcb6fb78 358void AliTRDseedV1::CookdEdx(Int_t nslices)
359{
360// Calculates average dE/dx for all slices and store them in the internal array fdEdx.
361//
362// Parameters:
363// nslices : number of slices for which dE/dx should be calculated
364// Output:
365// store results in the internal array fdEdx. This can be accessed with the method
366// AliTRDseedV1::GetdEdx()
367//
368// Detailed description
369// Calculates average dE/dx for all slices. Depending on the PID methode
370// the number of slices can be 3 (LQ) or 8(NN).
3ee48d6e 371// The calculation of dQ/dl are done using the tracklet fit results (see AliTRDseedV1::GetdQdl(Int_t))
bcb6fb78 372//
373// The following effects are included in the calculation:
374// 1. calibration values for t0 and vdrift (using x coordinate to calculate slice)
375// 2. cluster sharing (optional see AliTRDrecoParam::SetClusterSharing())
376// 3. cluster size
377//
378
e3cf3d02 379 Int_t nclusters[kNSlices];
380 memset(nclusters, 0, kNSlices*sizeof(Int_t));
381 memset(fdEdx, 0, kNSlices*sizeof(Float_t));
382
e73abf77 383 const Double_t kDriftLength = (.5 * AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
29b87567 384
3ee48d6e 385 AliTRDcluster *c = 0x0;
29b87567 386 for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++){
8e709c82 387 if(!(c = fClusters[ic]) && !(c = fClusters[ic+kNtb])) continue;
e73abf77 388 Float_t dx = TMath::Abs(fX0 - c->GetX());
29b87567 389
390 // Filter clusters for dE/dx calculation
391
392 // 1.consider calibration effects for slice determination
e73abf77 393 Int_t slice;
394 if(dx<kDriftLength){ // TODO should be replaced by c->IsInChamber()
395 slice = Int_t(dx * nslices / kDriftLength);
396 } else slice = c->GetX() < fX0 ? nslices-1 : 0;
397
398
29b87567 399 // 2. take sharing into account
3ee48d6e 400 Float_t w = c->IsShared() ? .5 : 1.;
29b87567 401
402 // 3. take into account large clusters TODO
403 //w *= c->GetNPads() > 3 ? .8 : 1.;
404
405 //CHECK !!!
406 fdEdx[slice] += w * GetdQdl(ic); //fdQdl[ic];
407 nclusters[slice]++;
408 } // End of loop over clusters
409
cd40b287 410 //if(fReconstructor->GetPIDMethod() == AliTRDReconstructor::kLQPID){
0d83b3a5 411 if(nslices == AliTRDpidUtil::kLQslices){
29b87567 412 // calculate mean charge per slice (only LQ PID)
413 for(int is=0; is<nslices; is++){
414 if(nclusters[is]) fdEdx[is] /= nclusters[is];
415 }
416 }
bcb6fb78 417}
418
e3cf3d02 419//_____________________________________________________________________________
420void AliTRDseedV1::CookLabels()
421{
422 //
423 // Cook 2 labels for seed
424 //
425
426 Int_t labels[200];
427 Int_t out[200];
428 Int_t nlab = 0;
429 for (Int_t i = 0; i < kNTimeBins; i++) {
430 if (!fClusters[i]) continue;
431 for (Int_t ilab = 0; ilab < 3; ilab++) {
432 if (fClusters[i]->GetLabel(ilab) >= 0) {
433 labels[nlab] = fClusters[i]->GetLabel(ilab);
434 nlab++;
435 }
436 }
437 }
438
fac58f00 439 fLabels[2] = AliMathBase::Freq(nlab,labels,out,kTRUE);
e3cf3d02 440 fLabels[0] = out[0];
441 if ((fLabels[2] > 1) && (out[3] > 1)) fLabels[1] = out[2];
442}
443
444
d937ad7a 445//____________________________________________________________________
446void AliTRDseedV1::GetClusterXY(const AliTRDcluster *c, Double_t &x, Double_t &y)
447{
448// Return corrected position of the cluster taking into
449// account variation of the drift velocity with drift length.
450
451
452 // drift velocity correction TODO to be moved to the clusterizer
453 const Float_t cx[] = {
454 -9.6280e-02, 1.3091e-01,-1.7415e-02,-9.9221e-02,-1.2040e-01,-9.5493e-02,
455 -5.0041e-02,-1.6726e-02, 3.5756e-03, 1.8611e-02, 2.6378e-02, 3.3823e-02,
456 3.4811e-02, 3.5282e-02, 3.5386e-02, 3.6047e-02, 3.5201e-02, 3.4384e-02,
457 3.2864e-02, 3.1932e-02, 3.2051e-02, 2.2539e-02,-2.5154e-02,-1.2050e-01,
458 -1.2050e-01
459 };
460
461 // PRF correction TODO to be replaced by the gaussian
462 // approximation with full error parametrization and // moved to the clusterizer
463 const Float_t cy[AliTRDgeometry::kNlayer][3] = {
464 { 4.014e-04, 8.605e-03, -6.880e+00},
465 {-3.061e-04, 9.663e-03, -6.789e+00},
466 { 1.124e-03, 1.105e-02, -6.825e+00},
467 {-1.527e-03, 1.231e-02, -6.777e+00},
468 { 2.150e-03, 1.387e-02, -6.783e+00},
469 {-1.296e-03, 1.486e-02, -6.825e+00}
470 };
471
472 Int_t ily = AliTRDgeometry::GetLayer(c->GetDetector());
473 x = c->GetX() - cx[c->GetLocalTimeBin()];
474 y = c->GetY() + cy[ily][0] + cy[ily][1] * TMath::Sin(cy[ily][2] * c->GetCenter());
475 return;
476}
b83573da 477
bcb6fb78 478//____________________________________________________________________
479Float_t AliTRDseedV1::GetdQdl(Int_t ic) const
480{
3ee48d6e 481// Using the linear approximation of the track inside one TRD chamber (TRD tracklet)
482// the charge per unit length can be written as:
483// BEGIN_LATEX
484// #frac{dq}{dl} = #frac{q_{c}}{dx * #sqrt{1 + #(){#frac{dy}{dx}}^{2}_{fit} + #(){#frac{dy}{dx}}^{2}_{ref}}}
485// END_LATEX
486// where qc is the total charge collected in the current time bin and dx is the length
487// of the time bin. For the moment (Jan 20 2009) only pad row cross corrections are
488// considered for the charge but none are applied for drift velocity variations along
489// the drift region or assymetry of the TRF
490//
491// Author : Alex Bercuci <A.Bercuci@gsi.de>
492//
493 Float_t dq = 0.;
494 if(fClusters[ic]) dq += TMath::Abs(fClusters[ic]->GetQ());
8e709c82 495 if(fClusters[ic+kNtb]) dq += TMath::Abs(fClusters[ic+kNtb]->GetQ());
496 if(dq<1.e-3 || fdX < 1.e-3) return 0.;
3ee48d6e 497
498 return dq/fdX/TMath::Sqrt(1. + fYfit[1]*fYfit[1] + fZref[1]*fZref[1]);
bcb6fb78 499}
500
0906e73e 501//____________________________________________________________________
e3cf3d02 502Float_t* AliTRDseedV1::GetProbability()
0906e73e 503{
504// Fill probability array for tracklet from the DB.
505//
506// Parameters
507//
508// Output
509// returns pointer to the probability array and 0x0 if missing DB access
510//
511// Detailed description
512
29b87567 513
514 // retrive calibration db
0906e73e 515 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
516 if (!calibration) {
517 AliError("No access to calibration data");
518 return 0x0;
519 }
520
3a039a31 521 if (!fReconstructor) {
522 AliError("Reconstructor not set.");
4ba1d6ae 523 return 0x0;
524 }
525
0906e73e 526 // Retrieve the CDB container class with the parametric detector response
3a039a31 527 const AliTRDCalPID *pd = calibration->GetPIDObject(fReconstructor->GetPIDMethod());
0906e73e 528 if (!pd) {
529 AliError("No access to AliTRDCalPID object");
530 return 0x0;
531 }
29b87567 532 //AliInfo(Form("Method[%d] : %s", fReconstructor->GetRecoParam() ->GetPIDMethod(), pd->IsA()->GetName()));
10f75631 533
29b87567 534 // calculate tracklet length TO DO
0906e73e 535 Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
536 /// TMath::Sqrt((1.0 - fSnp[iPlane]*fSnp[iPlane]) / (1.0 + fTgl[iPlane]*fTgl[iPlane]));
537
538 //calculate dE/dx
3a039a31 539 CookdEdx(fReconstructor->GetNdEdxSlices());
0906e73e 540
541 // Sets the a priori probabilities
542 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
ae4e8b84 543 fProb[ispec] = pd->GetProbability(ispec, fMom, &fdEdx[0], length, GetPlane());
0906e73e 544 }
545
29b87567 546 return &fProb[0];
0906e73e 547}
548
549//____________________________________________________________________
e4f2f73d 550Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
551{
552 //
553 // Returns a quality measurement of the current seed
554 //
555
e3cf3d02 556 Float_t zcorr = kZcorr ? fTilt * (fZfit[0] - fZref[0]) : 0.;
29b87567 557 return
558 .5 * TMath::Abs(18.0 - fN2)
559 + 10.* TMath::Abs(fYfit[1] - fYref[1])
560 + 5. * TMath::Abs(fYfit[0] - fYref[0] + zcorr)
e3cf3d02 561 + 2. * TMath::Abs(fZfit[0] - fZref[0]) / fPadLength;
e4f2f73d 562}
563
564//____________________________________________________________________
d937ad7a 565void AliTRDseedV1::GetCovAt(Double_t x, Double_t *cov) const
0906e73e 566{
d937ad7a 567// Computes covariance in the y-z plane at radial point x (in tracking coordinates)
568// and returns the results in the preallocated array cov[3] as :
569// cov[0] = Var(y)
570// cov[1] = Cov(yz)
571// cov[2] = Var(z)
572//
573// Details
574//
575// For the linear transformation
576// BEGIN_LATEX
577// Y = T_{x} X^{T}
578// END_LATEX
579// The error propagation has the general form
580// BEGIN_LATEX
581// C_{Y} = T_{x} C_{X} T_{x}^{T}
582// END_LATEX
583// We apply this formula 2 times. First to calculate the covariance of the tracklet
584// at point x we consider:
585// BEGIN_LATEX
586// T_{x} = (1 x); X=(y0 dy/dx); C_{X}=#(){#splitline{Var(y0) Cov(y0, dy/dx)}{Cov(y0, dy/dx) Var(dy/dx)}}
587// END_LATEX
588// and secondly to take into account the tilt angle
589// BEGIN_LATEX
590// T_{#alpha} = #(){#splitline{cos(#alpha) __ sin(#alpha)}{-sin(#alpha) __ cos(#alpha)}}; X=(y z); C_{X}=#(){#splitline{Var(y) 0}{0 Var(z)}}
591// END_LATEX
592//
593// using simple trigonometrics one can write for this last case
594// BEGIN_LATEX
595// 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})}}
596// END_LATEX
597// which can be aproximated for small alphas (2 deg) with
598// BEGIN_LATEX
599// 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}}}
600// END_LATEX
601//
602// before applying the tilt rotation we also apply systematic uncertainties to the tracklet
603// position which can be tunned from outside via the AliTRDrecoParam::SetSysCovMatrix(). They might
604// account for extra misalignment/miscalibration uncertainties.
605//
606// Author :
607// Alex Bercuci <A.Bercuci@gsi.de>
608// Date : Jan 8th 2009
609//
b1957d3c 610
611
d937ad7a 612 Double_t xr = fX0-x;
613 Double_t sy2 = fCov[0] +2.*xr*fCov[1] + xr*xr*fCov[2];
614 Double_t sz2 = fPadLength*fPadLength/12.;
0906e73e 615
d937ad7a 616 // insert systematic uncertainties
617 Double_t sys[15];
618 fReconstructor->GetRecoParam()->GetSysCovMatrix(sys);
619 sy2 += sys[0];
620 sz2 += sys[1];
621
622 // rotate covariance matrix
623 Double_t t2 = fTilt*fTilt;
624 Double_t correction = 1./(1. + t2);
625 cov[0] = (sy2+t2*sz2)*correction;
626 cov[1] = fTilt*(sz2 - sy2)*correction;
627 cov[2] = (t2*sy2+sz2)*correction;
628}
eb38ed55 629
0906e73e 630
d937ad7a 631//____________________________________________________________________
e3cf3d02 632void AliTRDseedV1::Calibrate()
d937ad7a 633{
e3cf3d02 634// Retrieve calibration and position parameters from OCDB.
635// The following information are used
d937ad7a 636// - detector index
e3cf3d02 637// - column and row position of first attached cluster. If no clusters are attached
638// to the tracklet a random central chamber position (c=70, r=7) will be used.
639//
640// The following information is cached in the tracklet
641// t0 (trigger delay)
642// drift velocity
643// PRF width
644// omega*tau = tg(a_L)
645// diffusion coefficients (longitudinal and transversal)
d937ad7a 646//
647// Author :
648// Alex Bercuci <A.Bercuci@gsi.de>
649// Date : Jan 8th 2009
650//
eb38ed55 651
d937ad7a 652 AliCDBManager *cdb = AliCDBManager::Instance();
653 if(cdb->GetRun() < 0){
654 AliError("OCDB manager not properly initialized");
655 return;
656 }
0906e73e 657
e3cf3d02 658 AliTRDcalibDB *calib = AliTRDcalibDB::Instance();
659 AliTRDCalROC *vdROC = calib->GetVdriftROC(fDet),
660 *t0ROC = calib->GetT0ROC(fDet);;
661 const AliTRDCalDet *vdDet = calib->GetVdriftDet();
662 const AliTRDCalDet *t0Det = calib->GetT0Det();
d937ad7a 663
664 Int_t col = 70, row = 7;
665 AliTRDcluster **c = &fClusters[0];
e3cf3d02 666 if(fN2){
d937ad7a 667 Int_t ic = 0;
e3cf3d02 668 while (ic<kNTimeBins && !(*c)){ic++; c++;}
d937ad7a 669 if(*c){
670 col = (*c)->GetPadCol();
671 row = (*c)->GetPadRow();
672 }
673 }
3a039a31 674
e3cf3d02 675 fT0 = t0Det->GetValue(fDet) + t0ROC->GetValue(col,row);
676 fVD = vdDet->GetValue(fDet) * vdROC->GetValue(col, row);
677 fS2PRF = calib->GetPRFWidth(fDet, col, row); fS2PRF *= fS2PRF;
678 fExB = AliTRDCommonParam::Instance()->GetOmegaTau(fVD);
679 AliTRDCommonParam::Instance()->GetDiffCoeff(fDiffL,
680 fDiffT, fVD);
681 SetBit(kCalib, kTRUE);
0906e73e 682}
683
0906e73e 684//____________________________________________________________________
29b87567 685void AliTRDseedV1::SetOwner()
0906e73e 686{
29b87567 687 //AliInfo(Form("own [%s] fOwner[%s]", own?"YES":"NO", fOwner?"YES":"NO"));
688
689 if(TestBit(kOwner)) return;
e3cf3d02 690 for(int ic=0; ic<kNTimeBins; ic++){
29b87567 691 if(!fClusters[ic]) continue;
692 fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
693 }
694 SetBit(kOwner);
0906e73e 695}
696
f29f13a6 697// //____________________________________________________________________
698// Bool_t AliTRDseedV1::AttachClustersIter(AliTRDtrackingChamber *chamber, Float_t quality, Bool_t kZcorr, AliTRDcluster *c)
699// {
700// //
701// // Iterative process to register clusters to the seed.
702// // In iteration 0 we try only one pad-row and if quality not
703// // sufficient we try 2 pad-rows (about 5% of tracks cross 2 pad-rows)
704// //
705// // debug level 7
706// //
707//
708// if(!fReconstructor->GetRecoParam() ){
709// AliError("Seed can not be used without a valid RecoParam.");
710// return kFALSE;
711// }
712//
713// AliTRDchamberTimeBin *layer = 0x0;
714// if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7){
715// AliTRDtrackingChamber ch(*chamber);
716// ch.SetOwner();
717// TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker);
718// cstreamer << "AttachClustersIter"
719// << "chamber.=" << &ch
720// << "tracklet.=" << this
721// << "\n";
722// }
723//
724// Float_t tquality;
725// Double_t kroady = fReconstructor->GetRecoParam() ->GetRoad1y();
726// Double_t kroadz = fPadLength * .5 + 1.;
727//
728// // initialize configuration parameters
729// Float_t zcorr = kZcorr ? fTilt * (fZfit[0] - fZref[0]) : 0.;
730// Int_t niter = kZcorr ? 1 : 2;
731//
732// Double_t yexp, zexp;
733// Int_t ncl = 0;
734// // start seed update
735// for (Int_t iter = 0; iter < niter; iter++) {
736// ncl = 0;
737// for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
738// if(!(layer = chamber->GetTB(iTime))) continue;
739// if(!Int_t(*layer)) continue;
740//
741// // define searching configuration
742// Double_t dxlayer = layer->GetX() - fX0;
743// if(c){
744// zexp = c->GetZ();
745// //Try 2 pad-rows in second iteration
746// if (iter > 0) {
747// zexp = fZref[0] + fZref[1] * dxlayer - zcorr;
748// if (zexp > c->GetZ()) zexp = c->GetZ() + fPadLength*0.5;
749// if (zexp < c->GetZ()) zexp = c->GetZ() - fPadLength*0.5;
750// }
751// } else zexp = fZref[0] + (kZcorr ? fZref[1] * dxlayer : 0.);
752// yexp = fYref[0] + fYref[1] * dxlayer - zcorr;
753//
754// // Get and register cluster
755// Int_t index = layer->SearchNearestCluster(yexp, zexp, kroady, kroadz);
756// if (index < 0) continue;
757// AliTRDcluster *cl = (*layer)[index];
758//
759// fIndexes[iTime] = layer->GetGlobalIndex(index);
760// fClusters[iTime] = cl;
761// // fY[iTime] = cl->GetY();
762// // fZ[iTime] = cl->GetZ();
763// ncl++;
764// }
765// if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d ncl [%d] = %d", iter, fDet, ncl));
766//
767// if(ncl>1){
768// // calculate length of the time bin (calibration aware)
769// Int_t irp = 0; Float_t x[2]={0., 0.}; Int_t tb[2] = {0,0};
e3cf3d02 770// for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
771// if(!fClusters[iTime]) continue;
f29f13a6 772// x[irp] = fClusters[iTime]->GetX();
773// tb[irp] = iTime;
774// irp++;
775// if(irp==2) break;
e3cf3d02 776// }
f29f13a6 777// Int_t dtb = tb[1] - tb[0];
778// fdX = dtb ? (x[0] - x[1]) / dtb : 0.15;
779//
780// // update X0 from the clusters (calibration/alignment aware)
781// for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
782// if(!(layer = chamber->GetTB(iTime))) continue;
783// if(!layer->IsT0()) continue;
784// if(fClusters[iTime]){
785// fX0 = fClusters[iTime]->GetX();
786// break;
787// } else { // we have to infere the position of the anode wire from the other clusters
788// for (Int_t jTime = iTime+1; jTime < AliTRDtrackerV1::GetNTimeBins(); jTime++) {
789// if(!fClusters[jTime]) continue;
790// fX0 = fClusters[jTime]->GetX() + fdX * (jTime - iTime);
791// break;
792// }
793// }
794// }
795//
796// // update YZ reference point
797// // TODO
798//
799// // update x reference positions (calibration/alignment aware)
800// // for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
801// // if(!fClusters[iTime]) continue;
802// // fX[iTime] = fX0 - fClusters[iTime]->GetX();
803// // }
804//
805// FitMI();
806// }
807// if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d nclFit [%d] = %d", iter, fDet, fN2));
808//
809// if(IsOK()){
810// tquality = GetQuality(kZcorr);
811// if(tquality < quality) break;
812// else quality = tquality;
813// }
814// kroadz *= 2.;
815// } // Loop: iter
816// if (!IsOK()) return kFALSE;
817//
818// if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=1) CookLabels();
819//
820// // load calibration params
821// Calibrate();
822// UpdateUsed();
823// return kTRUE;
824// }
e4f2f73d 825
826//____________________________________________________________________
b1957d3c 827Bool_t AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *chamber, Bool_t tilt)
e4f2f73d 828{
829 //
830 // Projective algorithm to attach clusters to seeding tracklets
831 //
832 // Parameters
833 //
834 // Output
835 //
836 // Detailed description
837 // 1. Collapse x coordinate for the full detector plane
838 // 2. truncated mean on y (r-phi) direction
839 // 3. purge clusters
840 // 4. truncated mean on z direction
841 // 5. purge clusters
842 // 6. fit tracklet
843 //
b1957d3c 844 Bool_t kPRINT = kFALSE;
29b87567 845 if(!fReconstructor->GetRecoParam() ){
846 AliError("Seed can not be used without a valid RecoParam.");
847 return kFALSE;
848 }
b1957d3c 849 // Initialize reco params for this tracklet
850 // 1. first time bin in the drift region
851 Int_t t0 = 4;
852 Int_t kClmin = Int_t(fReconstructor->GetRecoParam() ->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
29b87567 853
b1957d3c 854 Double_t syRef = TMath::Sqrt(fRefCov[0]);
29b87567 855 //define roads
b1957d3c 856 Double_t kroady = 1.;
857 //fReconstructor->GetRecoParam() ->GetRoad1y();
29b87567 858 Double_t kroadz = fPadLength * 1.5 + 1.;
b1957d3c 859 if(kPRINT) printf("AttachClusters() sy[%f] road[%f]\n", syRef, kroady);
29b87567 860
861 // working variables
b1957d3c 862 const Int_t kNrows = 16;
e3cf3d02 863 AliTRDcluster *clst[kNrows][kNTimeBins];
b1957d3c 864 Double_t cond[4], dx, dy, yt, zt,
e3cf3d02 865 yres[kNrows][kNTimeBins];
866 Int_t idxs[kNrows][kNTimeBins], ncl[kNrows], ncls = 0;
b1957d3c 867 memset(ncl, 0, kNrows*sizeof(Int_t));
e3cf3d02 868 memset(clst, 0, kNrows*kNTimeBins*sizeof(AliTRDcluster*));
b1957d3c 869
29b87567 870 // Do cluster projection
b1957d3c 871 AliTRDcluster *c = 0x0;
29b87567 872 AliTRDchamberTimeBin *layer = 0x0;
b1957d3c 873 Bool_t kBUFFER = kFALSE;
874 for (Int_t it = 0; it < AliTRDtrackerV1::GetNTimeBins(); it++) {
875 if(!(layer = chamber->GetTB(it))) continue;
29b87567 876 if(!Int_t(*layer)) continue;
877
b1957d3c 878 dx = fX0 - layer->GetX();
879 yt = fYref[0] - fYref[1] * dx;
880 zt = fZref[0] - fZref[1] * dx;
881 if(kPRINT) printf("\t%2d dx[%f] yt[%f] zt[%f]\n", it, dx, yt, zt);
882
883 // select clusters on a 5 sigmaKalman level
884 cond[0] = yt; cond[2] = kroady;
885 cond[1] = zt; cond[3] = kroadz;
886 Int_t n=0, idx[6];
887 layer->GetClusters(cond, idx, n, 6);
888 for(Int_t ic = n; ic--;){
889 c = (*layer)[idx[ic]];
890 dy = yt - c->GetY();
891 dy += tilt ? fTilt * (c->GetZ() - zt) : 0.;
892 // select clusters on a 3 sigmaKalman level
893/* if(tilt && TMath::Abs(dy) > 3.*syRef){
894 printf("too large !!!\n");
895 continue;
896 }*/
897 Int_t r = c->GetPadRow();
898 if(kPRINT) printf("\t\t%d dy[%f] yc[%f] r[%d]\n", ic, TMath::Abs(dy), c->GetY(), r);
899 clst[r][ncl[r]] = c;
900 idxs[r][ncl[r]] = idx[ic];
901 yres[r][ncl[r]] = dy;
902 ncl[r]++; ncls++;
903
e3cf3d02 904 if(ncl[r] >= kNTimeBins) {
905 AliWarning(Form("Cluster candidates reached limit %d. Some may be lost.", kNTimeBins));
b1957d3c 906 kBUFFER = kTRUE;
29b87567 907 break;
908 }
909 }
b1957d3c 910 if(kBUFFER) break;
29b87567 911 }
b1957d3c 912 if(kPRINT) printf("Found %d clusters\n", ncls);
913 if(ncls<kClmin) return kFALSE;
914
915 // analyze each row individualy
916 Double_t mean, syDis;
917 Int_t nrow[] = {0, 0, 0}, nr = 0, lr=-1;
918 for(Int_t ir=kNrows; ir--;){
919 if(!(ncl[ir])) continue;
920 if(lr>0 && lr-ir != 1){
921 if(kPRINT) printf("W - gap in rows attached !!\n");
29b87567 922 }
b1957d3c 923 if(kPRINT) printf("\tir[%d] lr[%d] n[%d]\n", ir, lr, ncl[ir]);
924 // Evaluate truncated mean on the y direction
925 if(ncl[ir] > 3) AliMathBase::EvaluateUni(ncl[ir], yres[ir], mean, syDis, Int_t(ncl[ir]*.8));
926 else {
927 mean = 0.; syDis = 0.;
928 }
929
930 // TODO check mean and sigma agains cluster resolution !!
931 if(kPRINT) printf("\tr[%2d] m[%f %5.3fsigma] s[%f]\n", ir, mean, TMath::Abs(mean/syRef), syDis);
932 // select clusters on a 3 sigmaDistr level
933 Bool_t kFOUND = kFALSE;
934 for(Int_t ic = ncl[ir]; ic--;){
935 if(yres[ir][ic] - mean > 3. * syDis){
936 clst[ir][ic] = 0x0; continue;
937 }
938 nrow[nr]++; kFOUND = kTRUE;
939 }
940 // exit loop
941 if(kFOUND) nr++;
942 lr = ir; if(nr>=3) break;
29b87567 943 }
b1957d3c 944 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]);
945
946 // classify cluster rows
947 Int_t row = -1;
948 switch(nr){
949 case 1:
950 row = lr;
951 break;
952 case 2:
953 SetBit(kRowCross, kTRUE); // mark pad row crossing
954 if(nrow[0] > nrow[1]){ row = lr+1; lr = -1;}
955 else{
956 row = lr; lr = 1;
957 nrow[2] = nrow[1];
958 nrow[1] = nrow[0];
959 nrow[0] = nrow[2];
29b87567 960 }
b1957d3c 961 break;
962 case 3:
963 SetBit(kRowCross, kTRUE); // mark pad row crossing
964 break;
29b87567 965 }
b1957d3c 966 if(kPRINT) printf("\trow[%d] n[%d]\n\n", row, nrow[0]);
967 if(row<0) return kFALSE;
29b87567 968
b1957d3c 969 // Select and store clusters
970 // We should consider here :
971 // 1. How far is the chamber boundary
972 // 2. How big is the mean
29b87567 973 fN2 = 0;
b1957d3c 974 for (Int_t ir = 0; ir < nr; ir++) {
975 Int_t jr = row + ir*lr;
976 if(kPRINT) printf("\tattach %d clusters for row %d\n", ncl[jr], jr);
977 for (Int_t ic = 0; ic < ncl[jr]; ic++) {
978 if(!(c = clst[jr][ic])) continue;
979 Int_t it = c->GetPadTime();
980 // TODO proper indexing of clusters !!
e3cf3d02 981 fIndexes[it+kNtb*ir] = chamber->GetTB(it)->GetGlobalIndex(idxs[jr][ic]);
982 fClusters[it+kNtb*ir] = c;
29b87567 983
b1957d3c 984 //printf("\tid[%2d] it[%d] idx[%d]\n", ic, it, fIndexes[it]);
985
986 fN2++;
987 }
988 }
989
29b87567 990 // number of minimum numbers of clusters expected for the tracklet
e4f2f73d 991 if (fN2 < kClmin){
29b87567 992 AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin));
e4f2f73d 993 fN2 = 0;
994 return kFALSE;
995 }
0906e73e 996
e3cf3d02 997 // Load calibration parameters for this tracklet
998 Calibrate();
b1957d3c 999
1000 // calculate dx for time bins in the drift region (calibration aware)
e3cf3d02 1001 Int_t irp = 0; Float_t x[2] = {0.,0.}; Int_t tb[2]={0,0};
b1957d3c 1002 for (Int_t it = t0; it < AliTRDtrackerV1::GetNTimeBins(); it++) {
1003 if(!fClusters[it]) continue;
1004 x[irp] = fClusters[it]->GetX();
1005 tb[irp] = it;
1006 irp++;
1007 if(irp==2) break;
e3cf3d02 1008 }
d86ed84c 1009 Int_t dtb = tb[1] - tb[0];
1010 fdX = dtb ? (x[0] - x[1]) / dtb : 0.15;
b1957d3c 1011
1012 // update X0 from the clusters (calibration/alignment aware) TODO remove dependence on x0 !!
1013 for (Int_t it = 0; it < AliTRDtrackerV1::GetNTimeBins(); it++) {
1014 if(!(layer = chamber->GetTB(it))) continue;
1015 if(!layer->IsT0()) continue;
1016 if(fClusters[it]){
1017 fX0 = fClusters[it]->GetX();
1018 break;
1019 } else { // we have to infere the position of the anode wire from the other clusters
1020 for (Int_t jt = it+1; jt < AliTRDtrackerV1::GetNTimeBins(); jt++) {
1021 if(!fClusters[jt]) continue;
1022 fX0 = fClusters[jt]->GetX() + fdX * (jt - it);
1023 break;
1024 }
1025 }
1026 }
1027
29b87567 1028 return kTRUE;
e4f2f73d 1029}
1030
03cef9b2 1031//____________________________________________________________
1032void AliTRDseedV1::Bootstrap(const AliTRDReconstructor *rec)
1033{
1034// Fill in all derived information. It has to be called after recovery from file or HLT.
1035// The primitive data are
1036// - list of clusters
1037// - detector (as the detector will be removed from clusters)
1038// - position of anode wire (fX0) - temporary
1039// - track reference position and direction
1040// - momentum of the track
1041// - time bin length [cm]
1042//
1043// A.Bercuci <A.Bercuci@gsi.de> Oct 30th 2008
1044//
1045 fReconstructor = rec;
1046 AliTRDgeometry g;
1047 AliTRDpadPlane *pp = g.GetPadPlane(fDet);
1048 fTilt = TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle());
1049 fPadLength = pp->GetLengthIPad();
e3cf3d02 1050 //fSnp = fYref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);
1051 //fTgl = fZref[1];
1052 fN2 = 0;// fMPads = 0.;
03cef9b2 1053 AliTRDcluster **cit = &fClusters[0];
e3cf3d02 1054 for(Int_t ic = kNTimeBins; ic--; cit++){
03cef9b2 1055 if(!(*cit)) return;
e3cf3d02 1056 fN2++;
1057/* fX[ic] = (*cit)->GetX() - fX0;
03cef9b2 1058 fY[ic] = (*cit)->GetY();
e3cf3d02 1059 fZ[ic] = (*cit)->GetZ();*/
03cef9b2 1060 }
e3cf3d02 1061 //Update(); //
1062 Fit();
03cef9b2 1063 CookLabels();
1064 GetProbability();
1065}
1066
1067
e4f2f73d 1068//____________________________________________________________________
d937ad7a 1069Bool_t AliTRDseedV1::Fit(Bool_t tilt, Int_t errors)
e4f2f73d 1070{
1071 //
1072 // Linear fit of the tracklet
1073 //
1074 // Parameters :
1075 //
1076 // Output :
1077 // True if successful
1078 //
1079 // Detailed description
1080 // 2. Check if tracklet crosses pad row boundary
1081 // 1. Calculate residuals in the y (r-phi) direction
1082 // 3. Do a Least Square Fit to the data
1083 //
1084
e3cf3d02 1085 if(!IsCalibrated()){
1086 AliWarning("Tracklet fit failed. Call Calibrate().");
1087 return kFALSE;
1088 }
1089
29b87567 1090 const Int_t kClmin = 8;
010d62b0 1091
9462866a 1092
1093 // cluster error parametrization parameters
010d62b0 1094 // 1. sy total charge
9462866a 1095 const Float_t sq0inv = 0.019962; // [1/q0]
1096 const Float_t sqb = 1.0281564; //[cm]
010d62b0 1097 // 2. sy for the PRF
1098 const Float_t scy[AliTRDgeometry::kNlayer][4] = {
d937ad7a 1099 {2.827e-02, 9.600e-04, 4.296e-01, 2.271e-02},
1100 {2.952e-02,-2.198e-04, 4.146e-01, 2.339e-02},
1101 {3.090e-02, 1.514e-03, 4.020e-01, 2.402e-02},
1102 {3.260e-02,-2.037e-03, 3.946e-01, 2.509e-02},
1103 {3.439e-02,-3.601e-04, 3.883e-01, 2.623e-02},
1104 {3.510e-02, 2.066e-03, 3.651e-01, 2.588e-02},
010d62b0 1105 };
1106 // 3. sy parallel to the track
d937ad7a 1107 const Float_t sy0 = 2.649e-02; // [cm]
1108 const Float_t sya = -8.864e-04; // [cm]
1109 const Float_t syb = -2.435e-01; // [cm]
1110
010d62b0 1111 // 4. sx parallel to the track
d937ad7a 1112 const Float_t sxgc = 5.427e-02;
1113 const Float_t sxgm = 7.783e-01;
1114 const Float_t sxgs = 2.743e-01;
1115 const Float_t sxe0 =-2.065e+00;
1116 const Float_t sxe1 =-2.978e-02;
1117
010d62b0 1118 // 5. sx perpendicular to the track
d937ad7a 1119// const Float_t sxd0 = 1.881e-02;
1120// const Float_t sxd1 =-4.101e-01;
1121// const Float_t sxd2 = 1.572e+00;
1122
2f7d6ac8 1123 // get track direction
1124 Double_t y0 = fYref[0];
1125 Double_t dydx = fYref[1];
1126 Double_t z0 = fZref[0];
1127 Double_t dzdx = fZref[1];
1128 Double_t yt, zt;
ae4e8b84 1129
29b87567 1130 const Int_t kNtb = AliTRDtrackerV1::GetNTimeBins();
e3cf3d02 1131 // calculation of tg^2(phi - a_L) and tg^2(a_L)
1132 Double_t tgg = (dydx-fExB)/(1.+dydx*fExB); tgg *= tgg;
1133 //Double_t exb2= fExB*fExB;
1134
b1957d3c 1135 //AliTRDtrackerV1::AliTRDLeastSquare fitterZ;
24d8660e 1136 TLinearFitter fitterY(1, "pol1");
29b87567 1137 // convertion factor from square to gauss distribution for sigma
b1957d3c 1138 //Double_t convert = 1./TMath::Sqrt(12.);
ae4e8b84 1139
29b87567 1140 // book cluster information
e3cf3d02 1141 Double_t qc[kNTimeBins], xc[kNTimeBins], yc[kNTimeBins], zc[kNTimeBins], sy[kNTimeBins];
1142
010d62b0 1143 Int_t ily = AliTRDgeometry::GetLayer(fDet);
e3cf3d02 1144 Int_t fN = 0;
9eb2d46c 1145 AliTRDcluster *c=0x0, **jc = &fClusters[0];
9eb2d46c 1146 for (Int_t ic=0; ic<kNtb; ic++, ++jc) {
b1957d3c 1147 //zRow[ic] = -1;
29b87567 1148 xc[ic] = -1.;
1149 yc[ic] = 999.;
1150 zc[ic] = 999.;
1151 sy[ic] = 0.;
b1957d3c 1152 //sz[ic] = 0.;
9eb2d46c 1153 if(!(c = (*jc))) continue;
29b87567 1154 if(!c->IsInChamber()) continue;
9462866a 1155
29b87567 1156 Float_t w = 1.;
1157 if(c->GetNPads()>4) w = .5;
1158 if(c->GetNPads()>5) w = .2;
010d62b0 1159
b1957d3c 1160 //zRow[fN] = c->GetPadRow();
e3cf3d02 1161 qc[fN] = TMath::Abs(c->GetQ());
d937ad7a 1162 // correct cluster position for PRF and v drift
e3cf3d02 1163 //Int_t jc = TMath::Max(fN-3, 0);
1164 //xc[fN] = c->GetXloc(fT0, fVD, &qc[jc], &xc[jc]/*, z0 - c->GetX()*dzdx*/);
1165 //Double_t s2 = fS2PRF + fDiffL*fDiffL*xc[fN]/(1.+2.*exb2)+tgg*xc[fN]*xc[fN]*exb2/12.;
1166 //yc[fN] = c->GetYloc(s2, fPadLength, xc[fN], fExB);
1167
1168 // uncalibrated cluster correction
1169 // TODO remove
d937ad7a 1170 Double_t x, y; GetClusterXY(c, x, y);
1171 xc[fN] = fX0 - x;
1172 yc[fN] = y;
2f7d6ac8 1173 zc[fN] = c->GetZ();
1174
1175 // extrapolated y value for the track
1176 yt = y0 - xc[fN]*dydx;
1177 // extrapolated z value for the track
1178 zt = z0 - xc[fN]*dzdx;
1179 // tilt correction
1180 if(tilt) yc[fN] -= fTilt*(zc[fN] - zt);
1181
010d62b0 1182 // ELABORATE CLUSTER ERROR
1183 // TODO to be moved to AliTRDcluster
010d62b0 1184 // basic y error (|| to track).
d937ad7a 1185 sy[fN] = xc[fN] < AliTRDgeometry::CamHght() ? 2. : sy0 + sya*TMath::Exp(1./(xc[fN]+syb));
1186 //printf("cluster[%d]\n\tsy[0] = %5.3e [um]\n", fN, sy[fN]*1.e4);
010d62b0 1187 // y error due to total charge
e3cf3d02 1188 sy[fN] += sqb*(1./qc[fN] - sq0inv);
d937ad7a 1189 //printf("\tsy[1] = %5.3e [um]\n", sy[fN]*1.e4);
010d62b0 1190 // y error due to PRF
1191 sy[fN] += scy[ily][0]*TMath::Gaus(c->GetCenter(), scy[ily][1], scy[ily][2]) - scy[ily][3];
d937ad7a 1192 //printf("\tsy[2] = %5.3e [um]\n", sy[fN]*1.e4);
1193
010d62b0 1194 sy[fN] *= sy[fN];
1195
1196 // ADD ERROR ON x
9462866a 1197 // error of drift length parallel to the track
d937ad7a 1198 Double_t sx = sxgc*TMath::Gaus(xc[fN], sxgm, sxgs) + TMath::Exp(sxe0+sxe1*xc[fN]); // [cm]
1199 //printf("\tsx[0] = %5.3e [um]\n", sx*1.e4);
9462866a 1200 // error of drift length perpendicular to the track
1201 //sx += sxd0 + sxd1*d + sxd2*d*d;
d937ad7a 1202 sx *= sx; // square sx
d937ad7a 1203
9462866a 1204 // add error from ExB
d937ad7a 1205 if(errors>0) sy[fN] += fExB*fExB*sx;
1206 //printf("\tsy[3] = %5.3e [um^2]\n", sy[fN]*1.e8);
1207
1208 // global radial error due to misalignment/miscalibration
1209 Double_t sx0 = 0.; sx0 *= sx0;
1210 // add sx contribution to sy due to track angle
1211 if(errors>1) sy[fN] += tgg*(sx+sx0);
1212 // TODO we should add tilt pad correction here
1213 //printf("\tsy[4] = %5.3e [um^2]\n", sy[fN]*1.e8);
1214 c->SetSigmaY2(sy[fN]);
1215
9462866a 1216 sy[fN] = TMath::Sqrt(sy[fN]);
e3cf3d02 1217 fitterY.AddPoint(&xc[fN], yc[fN], sy[fN]);
2f7d6ac8 1218 fN++;
29b87567 1219 }
47d5d320 1220 // to few clusters
2f7d6ac8 1221 if (fN < kClmin) return kFALSE;
1222
d937ad7a 1223 // fit XY
2f7d6ac8 1224 fitterY.Eval();
d937ad7a 1225 fYfit[0] = fitterY.GetParameter(0);
1226 fYfit[1] = -fitterY.GetParameter(1);
1227 // store covariance
1228 Double_t *p = fitterY.GetCovarianceMatrix();
1229 fCov[0] = p[0]; // variance of y0
1230 fCov[1] = p[1]; // covariance of y0, dydx
1231 fCov[2] = p[3]; // variance of dydx
b1957d3c 1232 // the ref radial position is set at the minimum of
1233 // the y variance of the tracklet
f29f13a6 1234 fX = -fCov[1]/fCov[2]; //fXref = fX0 - fXref;
1235 fS2Y = fCov[0] +2.*fX*fCov[1] + fX*fX*fCov[2];
b1957d3c 1236
1237 // fit XZ
1238 if(IsRowCross()){
1239 // TODO pad row cross position estimation !!!
1240 //AliInfo(Form("Padrow cross in detector %d", fDet));
1241 fZfit[0] = .5*(zc[0]+zc[fN-1]); fZfit[1] = 0.;
f29f13a6 1242 fS2Z = 0.02+1.55*fZref[1]; fS2Z *= fS2Z;
b1957d3c 1243 } else {
1244 fZfit[0] = zc[0]; fZfit[1] = 0.;
f29f13a6 1245 fS2Z = fPadLength*fPadLength/12.;
29b87567 1246 }
1247
29b87567 1248
b1957d3c 1249// // determine z offset of the fit
1250// Float_t zslope = 0.;
1251// Int_t nchanges = 0, nCross = 0;
1252// if(nz==2){ // tracklet is crossing pad row
1253// // Find the break time allowing one chage on pad-rows
1254// // with maximal number of accepted clusters
1255// Int_t padRef = zRow[0];
1256// for (Int_t ic=1; ic<fN; ic++) {
1257// if(zRow[ic] == padRef) continue;
1258//
1259// // debug
1260// if(zRow[ic-1] == zRow[ic]){
1261// printf("ERROR in pad row change!!!\n");
1262// }
1263//
1264// // evaluate parameters of the crossing point
1265// Float_t sx = (xc[ic-1] - xc[ic])*convert;
1266// fCross[0] = .5 * (xc[ic-1] + xc[ic]);
1267// fCross[2] = .5 * (zc[ic-1] + zc[ic]);
1268// fCross[3] = TMath::Max(dzdx * sx, .01);
1269// zslope = zc[ic-1] > zc[ic] ? 1. : -1.;
1270// padRef = zRow[ic];
1271// nCross = ic;
1272// nchanges++;
1273// }
1274// }
1275//
1276// // condition on nCross and reset nchanges TODO
1277//
1278// if(nchanges==1){
1279// if(dzdx * zslope < 0.){
1280// AliInfo("Tracklet-Track mismatch in dzdx. TODO.");
1281// }
1282//
1283//
1284// //zc[nc] = fitterZ.GetFunctionParameter(0);
1285// fCross[1] = fYfit[0] - fCross[0] * fYfit[1];
1286// fCross[0] = fX0 - fCross[0];
1287// }
29b87567 1288
29b87567 1289 return kTRUE;
e4f2f73d 1290}
1291
e4f2f73d 1292
f29f13a6 1293/*
e3cf3d02 1294//_____________________________________________________________________________
1295void AliTRDseedV1::FitMI()
1296{
1297//
1298// Fit the seed.
1299// Marian Ivanov's version
1300//
1301// linear fit on the y direction with respect to the reference direction.
1302// The residuals for each x (x = xc - x0) are deduced from:
1303// dy = y - yt (1)
1304// the tilting correction is written :
1305// y = yc + h*(zc-zt) (2)
1306// yt = y0+dy/dx*x (3)
1307// zt = z0+dz/dx*x (4)
1308// from (1),(2),(3) and (4)
1309// dy = yc - y0 - (dy/dx + h*dz/dx)*x + h*(zc-z0)
1310// the last term introduces the correction on y direction due to tilting pads. There are 2 ways to account for this:
1311// 1. use tilting correction for calculating the y
1312// 2. neglect tilting correction here and account for it in the error parametrization of the tracklet.
1313 const Float_t kRatio = 0.8;
1314 const Int_t kClmin = 5;
1315 const Float_t kmaxtan = 2;
1316
1317 if (TMath::Abs(fYref[1]) > kmaxtan){
1318 //printf("Exit: Abs(fYref[1]) = %3.3f, kmaxtan = %3.3f\n", TMath::Abs(fYref[1]), kmaxtan);
1319 return; // Track inclined too much
1320 }
1321
1322 Float_t sigmaexp = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction
1323 Float_t ycrosscor = fPadLength * fTilt * 0.5; // Y correction for crossing
1324 Int_t fNChange = 0;
1325
1326 Double_t sumw;
1327 Double_t sumwx;
1328 Double_t sumwx2;
1329 Double_t sumwy;
1330 Double_t sumwxy;
1331 Double_t sumwz;
1332 Double_t sumwxz;
1333
1334 // Buffering: Leave it constant fot Performance issues
1335 Int_t zints[kNtb]; // Histograming of the z coordinate
1336 // Get 1 and second max probable coodinates in z
1337 Int_t zouts[2*kNtb];
1338 Float_t allowedz[kNtb]; // Allowed z for given time bin
1339 Float_t yres[kNtb]; // Residuals from reference
1340 //Float_t anglecor = fTilt * fZref[1]; // Correction to the angle
1341
1342 Float_t pos[3*kNtb]; memset(pos, 0, 3*kNtb*sizeof(Float_t));
1343 Float_t *fX = &pos[0], *fY = &pos[kNtb], *fZ = &pos[2*kNtb];
1344
1345 Int_t fN = 0; AliTRDcluster *c = 0x0;
1346 fN2 = 0;
1347 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1348 yres[i] = 10000.0;
1349 if (!(c = fClusters[i])) continue;
1350 if(!c->IsInChamber()) continue;
1351 // Residual y
1352 //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] + fTilt*(fZ[i] - fZref[0]);
1353 fX[i] = fX0 - c->GetX();
1354 fY[i] = c->GetY();
1355 fZ[i] = c->GetZ();
1356 yres[i] = fY[i] - fTilt*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
1357 zints[fN] = Int_t(fZ[i]);
1358 fN++;
1359 }
1360
1361 if (fN < kClmin){
1362 //printf("Exit fN < kClmin: fN = %d\n", fN);
1363 return;
1364 }
1365 Int_t nz = AliTRDtrackerV1::Freq(fN, zints, zouts, kFALSE);
1366 Float_t fZProb = zouts[0];
1367 if (nz <= 1) zouts[3] = 0;
1368 if (zouts[1] + zouts[3] < kClmin) {
1369 //printf("Exit zouts[1] = %d, zouts[3] = %d\n",zouts[1],zouts[3]);
1370 return;
1371 }
1372
1373 // Z distance bigger than pad - length
1374 if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) zouts[3] = 0;
1375
1376 Int_t breaktime = -1;
1377 Bool_t mbefore = kFALSE;
1378 Int_t cumul[kNtb][2];
1379 Int_t counts[2] = { 0, 0 };
1380
1381 if (zouts[3] >= 3) {
1382
1383 //
1384 // Find the break time allowing one chage on pad-rows
1385 // with maximal number of accepted clusters
1386 //
1387 fNChange = 1;
1388 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1389 cumul[i][0] = counts[0];
1390 cumul[i][1] = counts[1];
1391 if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++;
1392 if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++;
1393 }
1394 Int_t maxcount = 0;
1395 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1396 Int_t after = cumul[AliTRDtrackerV1::GetNTimeBins()][0] - cumul[i][0];
1397 Int_t before = cumul[i][1];
1398 if (after + before > maxcount) {
1399 maxcount = after + before;
1400 breaktime = i;
1401 mbefore = kFALSE;
1402 }
1403 after = cumul[AliTRDtrackerV1::GetNTimeBins()-1][1] - cumul[i][1];
1404 before = cumul[i][0];
1405 if (after + before > maxcount) {
1406 maxcount = after + before;
1407 breaktime = i;
1408 mbefore = kTRUE;
1409 }
1410 }
1411 breaktime -= 1;
1412 }
1413
1414 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1415 if (i > breaktime) allowedz[i] = mbefore ? zouts[2] : zouts[0];
1416 if (i <= breaktime) allowedz[i] = (!mbefore) ? zouts[2] : zouts[0];
1417 }
1418
1419 if (((allowedz[0] > allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] < 0)) ||
1420 ((allowedz[0] < allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] > 0))) {
1421 //
1422 // Tracklet z-direction not in correspondance with track z direction
1423 //
1424 fNChange = 0;
1425 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1426 allowedz[i] = zouts[0]; // Only longest taken
1427 }
1428 }
1429
1430 if (fNChange > 0) {
1431 //
1432 // Cross pad -row tracklet - take the step change into account
1433 //
1434 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1435 if (!fClusters[i]) continue;
1436 if(!fClusters[i]->IsInChamber()) continue;
1437 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
1438 // Residual y
f29f13a6 1439 //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] + fTilt*(fZ[i] - fZref[0]);
e3cf3d02 1440 yres[i] = fY[i] - fTilt*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
f29f13a6 1441// if (TMath::Abs(fZ[i] - fZProb) > 2) {
1442// if (fZ[i] > fZProb) yres[i] += fTilt * fPadLength;
1443// if (fZ[i] < fZProb) yres[i] -= fTilt * fPadLength;
1444 }
e3cf3d02 1445 }
1446 }
1447
1448 Double_t yres2[kNtb];
1449 Double_t mean;
1450 Double_t sigma;
1451 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1452 if (!fClusters[i]) continue;
1453 if(!fClusters[i]->IsInChamber()) continue;
1454 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
1455 yres2[fN2] = yres[i];
1456 fN2++;
1457 }
1458 if (fN2 < kClmin) {
1459 //printf("Exit fN2 < kClmin: fN2 = %d\n", fN2);
1460 fN2 = 0;
1461 return;
1462 }
1463 AliMathBase::EvaluateUni(fN2,yres2,mean,sigma, Int_t(fN2*kRatio-2.));
1464 if (sigma < sigmaexp * 0.8) {
1465 sigma = sigmaexp;
1466 }
1467 //Float_t fSigmaY = sigma;
1468
1469 // Reset sums
1470 sumw = 0;
1471 sumwx = 0;
1472 sumwx2 = 0;
1473 sumwy = 0;
1474 sumwxy = 0;
1475 sumwz = 0;
1476 sumwxz = 0;
1477
1478 fN2 = 0;
1479 Float_t fMeanz = 0;
1480 Float_t fMPads = 0;
1481 fUsable = 0;
1482 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1483 if (!fClusters[i]) continue;
1484 if (!fClusters[i]->IsInChamber()) continue;
1485 if (TMath::Abs(fZ[i] - allowedz[i]) > 2){fClusters[i] = 0x0; continue;}
1486 if (TMath::Abs(yres[i] - mean) > 4.0 * sigma){fClusters[i] = 0x0; continue;}
1487 SETBIT(fUsable,i);
1488 fN2++;
1489 fMPads += fClusters[i]->GetNPads();
1490 Float_t weight = 1.0;
1491 if (fClusters[i]->GetNPads() > 4) weight = 0.5;
1492 if (fClusters[i]->GetNPads() > 5) weight = 0.2;
1493
1494
1495 Double_t x = fX[i];
1496 //printf("x = %7.3f dy = %7.3f fit %7.3f\n", x, yres[i], fY[i]-yres[i]);
1497
1498 sumw += weight;
1499 sumwx += x * weight;
1500 sumwx2 += x*x * weight;
1501 sumwy += weight * yres[i];
1502 sumwxy += weight * (yres[i]) * x;
1503 sumwz += weight * fZ[i];
1504 sumwxz += weight * fZ[i] * x;
1505
1506 }
1507
1508 if (fN2 < kClmin){
1509 //printf("Exit fN2 < kClmin(2): fN2 = %d\n",fN2);
1510 fN2 = 0;
1511 return;
1512 }
1513 fMeanz = sumwz / sumw;
1514 Float_t correction = 0;
1515 if (fNChange > 0) {
1516 // Tracklet on boundary
1517 if (fMeanz < fZProb) correction = ycrosscor;
1518 if (fMeanz > fZProb) correction = -ycrosscor;
1519 }
1520
1521 Double_t det = sumw * sumwx2 - sumwx * sumwx;
1522 fYfit[0] = (sumwx2 * sumwy - sumwx * sumwxy) / det;
1523 fYfit[1] = (sumw * sumwxy - sumwx * sumwy) / det;
1524
1525 fS2Y = 0;
1526 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1527 if (!TESTBIT(fUsable,i)) continue;
1528 Float_t delta = yres[i] - fYfit[0] - fYfit[1] * fX[i];
1529 fS2Y += delta*delta;
1530 }
1531 fS2Y = TMath::Sqrt(fS2Y / Float_t(fN2-2));
1532 // TEMPORARY UNTIL covariance properly calculated
1533 fS2Y = TMath::Max(fS2Y, Float_t(.1));
1534
1535 fZfit[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det;
1536 fZfit[1] = (sumw * sumwxz - sumwx * sumwz) / det;
1537// fYfitR[0] += fYref[0] + correction;
1538// fYfitR[1] += fYref[1];
1539// fYfit[0] = fYfitR[0];
1540 fYfit[1] = -fYfit[1];
1541
1542 UpdateUsed();
f29f13a6 1543}*/
e3cf3d02 1544
e4f2f73d 1545//___________________________________________________________________
203967fc 1546void AliTRDseedV1::Print(Option_t *o) const
e4f2f73d 1547{
1548 //
1549 // Printing the seedstatus
1550 //
1551
203967fc 1552 AliInfo(Form("Det[%3d] Tilt[%+6.2f] Pad[%5.2f]", fDet, fTilt, fPadLength));
f29f13a6 1553 AliInfo(Form("N[%2d] ", fN2));
203967fc 1554 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]));
1555 AliInfo(Form("Ref y[%7.2f] z[%7.2f] dydx[%5.2f] dzdx[%5.2f]", fYref[0], fZref[0], fYref[1], fZref[1]))
1556
1557
1558 if(strcmp(o, "a")!=0) return;
1559
4dc4dc2e 1560 AliTRDcluster* const* jc = &fClusters[0];
e3cf3d02 1561 for(int ic=0; ic<kNTimeBins; ic++, jc++) {
4dc4dc2e 1562 if(!(*jc)) continue;
203967fc 1563 (*jc)->Print(o);
4dc4dc2e 1564 }
e4f2f73d 1565}
47d5d320 1566
203967fc 1567
1568//___________________________________________________________________
1569Bool_t AliTRDseedV1::IsEqual(const TObject *o) const
1570{
1571 // Checks if current instance of the class has the same essential members
1572 // as the given one
1573
1574 if(!o) return kFALSE;
1575 const AliTRDseedV1 *inTracklet = dynamic_cast<const AliTRDseedV1*>(o);
1576 if(!inTracklet) return kFALSE;
1577
1578 for (Int_t i = 0; i < 2; i++){
e3cf3d02 1579 if ( fYref[i] != inTracklet->fYref[i] ) return kFALSE;
1580 if ( fZref[i] != inTracklet->fZref[i] ) return kFALSE;
203967fc 1581 }
1582
e3cf3d02 1583 if ( fS2Y != inTracklet->fS2Y ) return kFALSE;
1584 if ( fTilt != inTracklet->fTilt ) return kFALSE;
1585 if ( fPadLength != inTracklet->fPadLength ) return kFALSE;
203967fc 1586
e3cf3d02 1587 for (Int_t i = 0; i < kNTimeBins; i++){
1588// if ( fX[i] != inTracklet->GetX(i) ) return kFALSE;
1589// if ( fY[i] != inTracklet->GetY(i) ) return kFALSE;
1590// if ( fZ[i] != inTracklet->GetZ(i) ) return kFALSE;
1591 if ( fIndexes[i] != inTracklet->fIndexes[i] ) return kFALSE;
203967fc 1592 }
f29f13a6 1593// if ( fUsable != inTracklet->fUsable ) return kFALSE;
203967fc 1594
1595 for (Int_t i=0; i < 2; i++){
e3cf3d02 1596 if ( fYfit[i] != inTracklet->fYfit[i] ) return kFALSE;
1597 if ( fZfit[i] != inTracklet->fZfit[i] ) return kFALSE;
1598 if ( fLabels[i] != inTracklet->fLabels[i] ) return kFALSE;
203967fc 1599 }
1600
e3cf3d02 1601/* if ( fMeanz != inTracklet->GetMeanz() ) return kFALSE;
1602 if ( fZProb != inTracklet->GetZProb() ) return kFALSE;*/
1603 if ( fN2 != inTracklet->fN2 ) return kFALSE;
1604 if ( fNUsed != inTracklet->fNUsed ) return kFALSE;
1605 //if ( fFreq != inTracklet->GetFreq() ) return kFALSE;
1606 //if ( fNChange != inTracklet->GetNChange() ) return kFALSE;
203967fc 1607
e3cf3d02 1608 if ( fC != inTracklet->fC ) return kFALSE;
1609 //if ( fCC != inTracklet->GetCC() ) return kFALSE;
1610 if ( fChi2 != inTracklet->fChi2 ) return kFALSE;
203967fc 1611 // if ( fChi2Z != inTracklet->GetChi2Z() ) return kFALSE;
1612
e3cf3d02 1613 if ( fDet != inTracklet->fDet ) return kFALSE;
1614 if ( fMom != inTracklet->fMom ) return kFALSE;
1615 if ( fdX != inTracklet->fdX ) return kFALSE;
203967fc 1616
e3cf3d02 1617 for (Int_t iCluster = 0; iCluster < kNTimeBins; iCluster++){
203967fc 1618 AliTRDcluster *curCluster = fClusters[iCluster];
e3cf3d02 1619 AliTRDcluster *inCluster = inTracklet->fClusters[iCluster];
203967fc 1620 if (curCluster && inCluster){
1621 if (! curCluster->IsEqual(inCluster) ) {
1622 curCluster->Print();
1623 inCluster->Print();
1624 return kFALSE;
1625 }
1626 } else {
1627 // if one cluster exists, and corresponding
1628 // in other tracklet doesn't - return kFALSE
1629 if(curCluster || inCluster) return kFALSE;
1630 }
1631 }
1632 return kTRUE;
1633}