fix compiler warnings
[u/mrichter/AliRoot.git] / TRD / AliTRDseedV1.cxx
CommitLineData
e4f2f73d 1/**************************************************************************
29b87567 2* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3* *
4* Author: The ALICE Off-line Project. *
5* Contributors are mentioned in the code where appropriate. *
6* *
7* Permission to use, copy, modify and distribute this software and its *
8* documentation strictly for non-commercial purposes is hereby granted *
9* without fee, provided that the above copyright notice appears in all *
10* copies and that both the copyright notice and this permission notice *
11* appear in the supporting documentation. The authors make no claims *
12* about the suitability of this software for any purpose. It is *
13* provided "as is" without express or implied warranty. *
14**************************************************************************/
e4f2f73d 15
16/* $Id$ */
17
18////////////////////////////////////////////////////////////////////////////
dd8059a8 19////
20// The TRD offline tracklet
21//
22// The running horse of the TRD reconstruction. The following tasks are preformed:
23// 1. Clusters attachment to tracks based on prior information stored at tracklet level (see AttachClusters)
24// 2. Clusters position recalculation based on track information (see GetClusterXY and Fit)
25// 3. Cluster error parametrization recalculation (see Fit)
26// 4. Linear track approximation (Fit)
27// 5. Optimal position (including z estimate for pad row cross tracklets) and covariance matrix of the track fit inside one TRD chamber (Fit)
28// 6. Tilt pad correction and systematic effects (GetCovAt)
29// 7. dEdx calculation (CookdEdx)
30// 8. PID probabilities estimation (CookPID)
31//
e4f2f73d 32// Authors: //
33// Alex Bercuci <A.Bercuci@gsi.de> //
34// Markus Fasel <M.Fasel@gsi.de> //
35// //
36////////////////////////////////////////////////////////////////////////////
37
38#include "TMath.h"
39#include "TLinearFitter.h"
eb38ed55 40#include "TClonesArray.h" // tmp
41#include <TTreeStream.h>
e4f2f73d 42
43#include "AliLog.h"
44#include "AliMathBase.h"
d937ad7a 45#include "AliCDBManager.h"
46#include "AliTracker.h"
e4f2f73d 47
03cef9b2 48#include "AliTRDpadPlane.h"
e4f2f73d 49#include "AliTRDcluster.h"
f3d3af1b 50#include "AliTRDseedV1.h"
51#include "AliTRDtrackV1.h"
e4f2f73d 52#include "AliTRDcalibDB.h"
eb38ed55 53#include "AliTRDchamberTimeBin.h"
54#include "AliTRDtrackingChamber.h"
55#include "AliTRDtrackerV1.h"
e4f2f73d 56#include "AliTRDrecoParam.h"
a076fc2f 57#include "AliTRDCommonParam.h"
d937ad7a 58
0906e73e 59#include "Cal/AliTRDCalPID.h"
d937ad7a 60#include "Cal/AliTRDCalROC.h"
61#include "Cal/AliTRDCalDet.h"
e4f2f73d 62
e4f2f73d 63ClassImp(AliTRDseedV1)
64
4d6aee34 65TLinearFitter *AliTRDseedV1::fgFitterY = NULL;
66TLinearFitter *AliTRDseedV1::fgFitterZ = NULL;
f301a656 67
e4f2f73d 68//____________________________________________________________________
ae4e8b84 69AliTRDseedV1::AliTRDseedV1(Int_t det)
3e778975 70 :AliTRDtrackletBase()
4d6aee34 71 ,fkReconstructor(NULL)
72 ,fClusterIter(NULL)
e3cf3d02 73 ,fExB(0.)
74 ,fVD(0.)
75 ,fT0(0.)
76 ,fS2PRF(0.)
77 ,fDiffL(0.)
78 ,fDiffT(0.)
ae4e8b84 79 ,fClusterIdx(0)
7c3eecb8 80 ,fErrorMsg(0)
3e778975 81 ,fN(0)
ae4e8b84 82 ,fDet(det)
b25a5e9e 83 ,fPt(0.)
bcb6fb78 84 ,fdX(0.)
e3cf3d02 85 ,fX0(0.)
86 ,fX(0.)
87 ,fY(0.)
88 ,fZ(0.)
89 ,fS2Y(0.)
90 ,fS2Z(0.)
91 ,fC(0.)
92 ,fChi2(0.)
e4f2f73d 93{
94 //
95 // Constructor
96 //
f301a656 97 memset(fIndexes,0xFF,kNclusters*sizeof(fIndexes[0]));
8d2bec9e 98 memset(fClusters, 0, kNclusters*sizeof(AliTRDcluster*));
dd8059a8 99 memset(fPad, 0, 3*sizeof(Float_t));
e3cf3d02 100 fYref[0] = 0.; fYref[1] = 0.;
101 fZref[0] = 0.; fZref[1] = 0.;
102 fYfit[0] = 0.; fYfit[1] = 0.;
103 fZfit[0] = 0.; fZfit[1] = 0.;
8d2bec9e 104 memset(fdEdx, 0, kNslices*sizeof(Float_t));
29b87567 105 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = -1.;
e3cf3d02 106 fLabels[0]=-1; fLabels[1]=-1; // most freq MC labels
107 fLabels[2]=0; // number of different labels for tracklet
16cca13f 108 memset(fRefCov, 0, 7*sizeof(Double_t));
d937ad7a 109 // covariance matrix [diagonal]
110 // default sy = 200um and sz = 2.3 cm
111 fCov[0] = 4.e-4; fCov[1] = 0.; fCov[2] = 5.3;
f29f13a6 112 SetStandAlone(kFALSE);
e4f2f73d 113}
114
115//____________________________________________________________________
0906e73e 116AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
3e778975 117 :AliTRDtrackletBase((AliTRDtrackletBase&)ref)
4d6aee34 118 ,fkReconstructor(NULL)
119 ,fClusterIter(NULL)
e3cf3d02 120 ,fExB(0.)
121 ,fVD(0.)
122 ,fT0(0.)
123 ,fS2PRF(0.)
124 ,fDiffL(0.)
125 ,fDiffT(0.)
ae4e8b84 126 ,fClusterIdx(0)
7c3eecb8 127 ,fErrorMsg(0)
3e778975 128 ,fN(0)
e3cf3d02 129 ,fDet(-1)
b25a5e9e 130 ,fPt(0.)
e3cf3d02 131 ,fdX(0.)
132 ,fX0(0.)
133 ,fX(0.)
134 ,fY(0.)
135 ,fZ(0.)
136 ,fS2Y(0.)
137 ,fS2Z(0.)
138 ,fC(0.)
139 ,fChi2(0.)
e4f2f73d 140{
141 //
142 // Copy Constructor performing a deep copy
143 //
e3cf3d02 144 if(this != &ref){
145 ref.Copy(*this);
146 }
29b87567 147 SetBit(kOwner, kFALSE);
f29f13a6 148 SetStandAlone(ref.IsStandAlone());
fbb2ea06 149}
d9950a5a 150
0906e73e 151
e4f2f73d 152//____________________________________________________________________
153AliTRDseedV1& AliTRDseedV1::operator=(const AliTRDseedV1 &ref)
154{
155 //
156 // Assignment Operator using the copy function
157 //
158
29b87567 159 if(this != &ref){
160 ref.Copy(*this);
161 }
221ab7e0 162 SetBit(kOwner, kFALSE);
163
29b87567 164 return *this;
e4f2f73d 165}
166
167//____________________________________________________________________
168AliTRDseedV1::~AliTRDseedV1()
169{
170 //
171 // Destructor. The RecoParam object belongs to the underlying tracker.
172 //
173
29b87567 174 //printf("I-AliTRDseedV1::~AliTRDseedV1() : Owner[%s]\n", IsOwner()?"YES":"NO");
e4f2f73d 175
e3cf3d02 176 if(IsOwner()) {
8d2bec9e 177 for(int itb=0; itb<kNclusters; itb++){
29b87567 178 if(!fClusters[itb]) continue;
179 //AliInfo(Form("deleting c %p @ %d", fClusters[itb], itb));
180 delete fClusters[itb];
4d6aee34 181 fClusters[itb] = NULL;
29b87567 182 }
e3cf3d02 183 }
e4f2f73d 184}
185
186//____________________________________________________________________
187void AliTRDseedV1::Copy(TObject &ref) const
188{
189 //
190 // Copy function
191 //
192
29b87567 193 //AliInfo("");
194 AliTRDseedV1 &target = (AliTRDseedV1 &)ref;
195
4d6aee34 196 target.fkReconstructor = fkReconstructor;
197 target.fClusterIter = NULL;
e3cf3d02 198 target.fExB = fExB;
199 target.fVD = fVD;
200 target.fT0 = fT0;
201 target.fS2PRF = fS2PRF;
202 target.fDiffL = fDiffL;
203 target.fDiffT = fDiffT;
ae4e8b84 204 target.fClusterIdx = 0;
7c3eecb8 205 target.fErrorMsg = fErrorMsg;
3e778975 206 target.fN = fN;
ae4e8b84 207 target.fDet = fDet;
b25a5e9e 208 target.fPt = fPt;
29b87567 209 target.fdX = fdX;
e3cf3d02 210 target.fX0 = fX0;
211 target.fX = fX;
212 target.fY = fY;
213 target.fZ = fZ;
214 target.fS2Y = fS2Y;
215 target.fS2Z = fS2Z;
216 target.fC = fC;
217 target.fChi2 = fChi2;
29b87567 218
8d2bec9e 219 memcpy(target.fIndexes, fIndexes, kNclusters*sizeof(Int_t));
220 memcpy(target.fClusters, fClusters, kNclusters*sizeof(AliTRDcluster*));
dd8059a8 221 memcpy(target.fPad, fPad, 3*sizeof(Float_t));
e3cf3d02 222 target.fYref[0] = fYref[0]; target.fYref[1] = fYref[1];
223 target.fZref[0] = fZref[0]; target.fZref[1] = fZref[1];
224 target.fYfit[0] = fYfit[0]; target.fYfit[1] = fYfit[1];
225 target.fZfit[0] = fZfit[0]; target.fZfit[1] = fZfit[1];
8d2bec9e 226 memcpy(target.fdEdx, fdEdx, kNslices*sizeof(Float_t));
e3cf3d02 227 memcpy(target.fProb, fProb, AliPID::kSPECIES*sizeof(Float_t));
228 memcpy(target.fLabels, fLabels, 3*sizeof(Int_t));
16cca13f 229 memcpy(target.fRefCov, fRefCov, 7*sizeof(Double_t));
e3cf3d02 230 memcpy(target.fCov, fCov, 3*sizeof(Double_t));
29b87567 231
e3cf3d02 232 TObject::Copy(ref);
e4f2f73d 233}
234
0906e73e 235
236//____________________________________________________________
f3d3af1b 237Bool_t AliTRDseedV1::Init(AliTRDtrackV1 *track)
0906e73e 238{
239// Initialize this tracklet using the track information
240//
241// Parameters:
242// track - the TRD track used to initialize the tracklet
243//
244// Detailed description
245// The function sets the starting point and direction of the
246// tracklet according to the information from the TRD track.
247//
248// Caution
249// The TRD track has to be propagated to the beginning of the
250// chamber where the tracklet will be constructed
251//
252
29b87567 253 Double_t y, z;
254 if(!track->GetProlongation(fX0, y, z)) return kFALSE;
16cca13f 255 Update(track);
29b87567 256 return kTRUE;
0906e73e 257}
258
bcb6fb78 259
e3cf3d02 260//_____________________________________________________________________________
261void AliTRDseedV1::Reset()
262{
263 //
264 // Reset seed
265 //
266 fExB=0.;fVD=0.;fT0=0.;fS2PRF=0.;
267 fDiffL=0.;fDiffT=0.;
3e778975 268 fClusterIdx=0;
7c3eecb8 269 fErrorMsg = 0;
3e778975 270 fN=0;
dd8059a8 271 fDet=-1;
b25a5e9e 272 fPt=0.;
e3cf3d02 273 fdX=0.;fX0=0.; fX=0.; fY=0.; fZ=0.;
274 fS2Y=0.; fS2Z=0.;
275 fC=0.; fChi2 = 0.;
276
8d2bec9e 277 for(Int_t ic=kNclusters; ic--;) fIndexes[ic] = -1;
278 memset(fClusters, 0, kNclusters*sizeof(AliTRDcluster*));
dd8059a8 279 memset(fPad, 0, 3*sizeof(Float_t));
e3cf3d02 280 fYref[0] = 0.; fYref[1] = 0.;
281 fZref[0] = 0.; fZref[1] = 0.;
282 fYfit[0] = 0.; fYfit[1] = 0.;
283 fZfit[0] = 0.; fZfit[1] = 0.;
8d2bec9e 284 memset(fdEdx, 0, kNslices*sizeof(Float_t));
e3cf3d02 285 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = -1.;
286 fLabels[0]=-1; fLabels[1]=-1; // most freq MC labels
287 fLabels[2]=0; // number of different labels for tracklet
16cca13f 288 memset(fRefCov, 0, 7*sizeof(Double_t));
e3cf3d02 289 // covariance matrix [diagonal]
290 // default sy = 200um and sz = 2.3 cm
291 fCov[0] = 4.e-4; fCov[1] = 0.; fCov[2] = 5.3;
292}
293
bcb6fb78 294//____________________________________________________________________
16cca13f 295void AliTRDseedV1::Update(const AliTRDtrackV1 *trk)
b1957d3c 296{
297 // update tracklet reference position from the TRD track
b1957d3c 298
e3cf3d02 299 Double_t fSnp = trk->GetSnp();
300 Double_t fTgl = trk->GetTgl();
b25a5e9e 301 fPt = trk->Pt();
1fd9389f 302 Double_t norm =1./TMath::Sqrt(1. - fSnp*fSnp);
303 fYref[1] = fSnp*norm;
304 fZref[1] = fTgl*norm;
b1957d3c 305 SetCovRef(trk->GetCovariance());
306
307 Double_t dx = trk->GetX() - fX0;
308 fYref[0] = trk->GetY() - dx*fYref[1];
309 fZref[0] = trk->GetZ() - dx*fZref[1];
310}
311
e3cf3d02 312//_____________________________________________________________________________
313void AliTRDseedV1::UpdateUsed()
314{
315 //
f29f13a6 316 // Calculate number of used clusers in the tracklet
e3cf3d02 317 //
318
3e778975 319 Int_t nused = 0, nshared = 0;
8d2bec9e 320 for (Int_t i = kNclusters; i--; ) {
e3cf3d02 321 if (!fClusters[i]) continue;
3e778975 322 if(fClusters[i]->IsUsed()){
323 nused++;
324 } else if(fClusters[i]->IsShared()){
325 if(IsStandAlone()) nused++;
326 else nshared++;
327 }
e3cf3d02 328 }
3e778975 329 SetNUsed(nused);
330 SetNShared(nshared);
e3cf3d02 331}
332
333//_____________________________________________________________________________
334void AliTRDseedV1::UseClusters()
335{
336 //
337 // Use clusters
338 //
f29f13a6 339 // In stand alone mode:
340 // Clusters which are marked as used or shared from another track are
341 // removed from the tracklet
342 //
343 // In barrel mode:
344 // - Clusters which are used by another track become shared
345 // - Clusters which are attached to a kink track become shared
346 //
e3cf3d02 347 AliTRDcluster **c = &fClusters[0];
8d2bec9e 348 for (Int_t ic=kNclusters; ic--; c++) {
e3cf3d02 349 if(!(*c)) continue;
f29f13a6 350 if(IsStandAlone()){
351 if((*c)->IsShared() || (*c)->IsUsed()){
b82b4de1 352 if((*c)->IsShared()) SetNShared(GetNShared()-1);
353 else SetNUsed(GetNUsed()-1);
4d6aee34 354 (*c) = NULL;
f29f13a6 355 fIndexes[ic] = -1;
3e778975 356 SetN(GetN()-1);
3e778975 357 continue;
f29f13a6 358 }
3e778975 359 } else {
f29f13a6 360 if((*c)->IsUsed() || IsKink()){
3e778975 361 (*c)->SetShared();
362 continue;
f29f13a6 363 }
364 }
365 (*c)->Use();
e3cf3d02 366 }
367}
368
369
f29f13a6 370
b1957d3c 371//____________________________________________________________________
bcb6fb78 372void AliTRDseedV1::CookdEdx(Int_t nslices)
373{
374// Calculates average dE/dx for all slices and store them in the internal array fdEdx.
375//
376// Parameters:
377// nslices : number of slices for which dE/dx should be calculated
378// Output:
379// store results in the internal array fdEdx. This can be accessed with the method
380// AliTRDseedV1::GetdEdx()
381//
382// Detailed description
383// Calculates average dE/dx for all slices. Depending on the PID methode
384// the number of slices can be 3 (LQ) or 8(NN).
3ee48d6e 385// The calculation of dQ/dl are done using the tracklet fit results (see AliTRDseedV1::GetdQdl(Int_t))
bcb6fb78 386//
387// The following effects are included in the calculation:
388// 1. calibration values for t0 and vdrift (using x coordinate to calculate slice)
389// 2. cluster sharing (optional see AliTRDrecoParam::SetClusterSharing())
390// 3. cluster size
391//
392
8d2bec9e 393 Int_t nclusters[kNslices];
394 memset(nclusters, 0, kNslices*sizeof(Int_t));
395 memset(fdEdx, 0, kNslices*sizeof(Float_t));
e3cf3d02 396
e73abf77 397 const Double_t kDriftLength = (.5 * AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
29b87567 398
4d6aee34 399 AliTRDcluster *c = NULL;
29b87567 400 for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++){
8e709c82 401 if(!(c = fClusters[ic]) && !(c = fClusters[ic+kNtb])) continue;
e73abf77 402 Float_t dx = TMath::Abs(fX0 - c->GetX());
29b87567 403
404 // Filter clusters for dE/dx calculation
405
406 // 1.consider calibration effects for slice determination
e73abf77 407 Int_t slice;
408 if(dx<kDriftLength){ // TODO should be replaced by c->IsInChamber()
409 slice = Int_t(dx * nslices / kDriftLength);
410 } else slice = c->GetX() < fX0 ? nslices-1 : 0;
411
412
29b87567 413 // 2. take sharing into account
3e778975 414 Float_t w = /*c->IsShared() ? .5 :*/ 1.;
29b87567 415
416 // 3. take into account large clusters TODO
417 //w *= c->GetNPads() > 3 ? .8 : 1.;
418
419 //CHECK !!!
420 fdEdx[slice] += w * GetdQdl(ic); //fdQdl[ic];
421 nclusters[slice]++;
422 } // End of loop over clusters
423
4d6aee34 424 //if(fkReconstructor->GetPIDMethod() == AliTRDReconstructor::kLQPID){
0d83b3a5 425 if(nslices == AliTRDpidUtil::kLQslices){
29b87567 426 // calculate mean charge per slice (only LQ PID)
427 for(int is=0; is<nslices; is++){
428 if(nclusters[is]) fdEdx[is] /= nclusters[is];
429 }
430 }
bcb6fb78 431}
432
e3cf3d02 433//_____________________________________________________________________________
434void AliTRDseedV1::CookLabels()
435{
436 //
437 // Cook 2 labels for seed
438 //
439
440 Int_t labels[200];
441 Int_t out[200];
442 Int_t nlab = 0;
8d2bec9e 443 for (Int_t i = 0; i < kNclusters; i++) {
e3cf3d02 444 if (!fClusters[i]) continue;
445 for (Int_t ilab = 0; ilab < 3; ilab++) {
446 if (fClusters[i]->GetLabel(ilab) >= 0) {
447 labels[nlab] = fClusters[i]->GetLabel(ilab);
448 nlab++;
449 }
450 }
451 }
452
fac58f00 453 fLabels[2] = AliMathBase::Freq(nlab,labels,out,kTRUE);
e3cf3d02 454 fLabels[0] = out[0];
455 if ((fLabels[2] > 1) && (out[3] > 1)) fLabels[1] = out[2];
456}
457
458
d937ad7a 459//____________________________________________________________________
0b433f72 460Float_t AliTRDseedV1::GetdQdl(Int_t ic, Float_t *dl) const
bcb6fb78 461{
3ee48d6e 462// Using the linear approximation of the track inside one TRD chamber (TRD tracklet)
463// the charge per unit length can be written as:
464// BEGIN_LATEX
500851ab 465// #frac{dq}{dl} = #frac{q_{c}}{dx * #sqrt{1 + #(){#frac{dy}{dx}}^{2}_{fit} + #(){#frac{dz}{dx}}^{2}_{ref}}}
3ee48d6e 466// END_LATEX
467// where qc is the total charge collected in the current time bin and dx is the length
0b433f72 468// of the time bin.
469// The following correction are applied :
470// - charge : pad row cross corrections
471// [diffusion and TRF assymetry] TODO
472// - dx : anisochronity, track inclination - see Fit and AliTRDcluster::GetXloc()
473// and AliTRDcluster::GetYloc() for the effects taken into account
3ee48d6e 474//
0fa1a8ee 475//Begin_Html
476//<img src="TRD/trackletDQDT.gif">
477//End_Html
478// In the picture the energy loss measured on the tracklet as a function of drift time [left] and respectively
479// drift length [right] for different particle species is displayed.
3ee48d6e 480// Author : Alex Bercuci <A.Bercuci@gsi.de>
481//
482 Float_t dq = 0.;
5d401b45 483 // check whether both clusters are inside the chamber
484 Bool_t hasClusterInChamber = kFALSE;
485 if(fClusters[ic] && fClusters[ic]->IsInChamber()){
486 hasClusterInChamber = kTRUE;
1742f24c 487 dq += TMath::Abs(fClusters[ic]->GetQ());
5d401b45 488 }else if(fClusters[ic+kNtb] && fClusters[ic+kNtb]->IsInChamber()){
489 hasClusterInChamber = kTRUE;
490 dq += TMath::Abs(fClusters[ic+kNtb]->GetQ());
1742f24c 491 }
5d401b45 492 if(!hasClusterInChamber) return 0.;
0b433f72 493 if(dq<1.e-3) return 0.;
3ee48d6e 494
a2abcbc5 495 Double_t dx = fdX;
496 if(ic-1>=0 && ic+1<kNtb){
497 Float_t x2(0.), x1(0.);
5d401b45 498 // try to estimate upper radial position (find the cluster which is inside the chamber)
499 if(fClusters[ic-1] && fClusters[ic-1]->IsInChamber()) x2 = fClusters[ic-1]->GetX();
500 else if(fClusters[ic-1+kNtb] && fClusters[ic-1+kNtb]->IsInChamber()) x2 = fClusters[ic-1+kNtb]->GetX();
501 else if(fClusters[ic] && fClusters[ic]->IsInChamber()) x2 = fClusters[ic]->GetX()+fdX;
a2abcbc5 502 else x2 = fClusters[ic+kNtb]->GetX()+fdX;
5d401b45 503 // try to estimate lower radial position (find the cluster which is inside the chamber)
504 if(fClusters[ic+1] && fClusters[ic+1]->IsInChamber()) x1 = fClusters[ic+1]->GetX();
505 else if(fClusters[ic+1+kNtb] && fClusters[ic+1+kNtb]->IsInChamber()) x1 = fClusters[ic+1+kNtb]->GetX();
506 else if(fClusters[ic] && fClusters[ic]->IsInChamber()) x1 = fClusters[ic]->GetX()-fdX;
a2abcbc5 507 else x1 = fClusters[ic+kNtb]->GetX()-fdX;
508
509 dx = .5*(x2 - x1);
510 }
0b433f72 511 dx *= TMath::Sqrt(1. + fYfit[1]*fYfit[1] + fZref[1]*fZref[1]);
0b433f72 512 if(dl) (*dl) = dx;
513 return dq/dx;
bcb6fb78 514}
515
0b433f72 516//____________________________________________________________
517Float_t AliTRDseedV1::GetMomentum(Float_t *err) const
518{
519// Returns momentum of the track after update with the current tracklet as:
520// BEGIN_LATEX
521// p=#frac{1}{1/p_{t}} #sqrt{1+tgl^{2}}
522// END_LATEX
523// and optionally the momentum error (if err is not null).
524// The estimated variance of the momentum is given by:
525// BEGIN_LATEX
526// #sigma_{p}^{2} = (#frac{dp}{dp_{t}})^{2} #sigma_{p_{t}}^{2}+(#frac{dp}{dtgl})^{2} #sigma_{tgl}^{2}+2#frac{dp}{dp_{t}}#frac{dp}{dtgl} cov(tgl,1/p_{t})
527// END_LATEX
528// which can be simplified to
529// BEGIN_LATEX
530// #sigma_{p}^{2} = p^{2}p_{t}^{4}tgl^{2}#sigma_{tgl}^{2}-2p^{2}p_{t}^{3}tgl cov(tgl,1/p_{t})+p^{2}p_{t}^{2}#sigma_{1/p_{t}}^{2}
531// END_LATEX
532//
533
534 Double_t p = fPt*TMath::Sqrt(1.+fZref[1]*fZref[1]);
535 Double_t p2 = p*p;
536 Double_t tgl2 = fZref[1]*fZref[1];
537 Double_t pt2 = fPt*fPt;
538 if(err){
539 Double_t s2 =
540 p2*tgl2*pt2*pt2*fRefCov[4]
541 -2.*p2*fZref[1]*fPt*pt2*fRefCov[5]
542 +p2*pt2*fRefCov[6];
543 (*err) = TMath::Sqrt(s2);
544 }
545 return p;
546}
547
548
0906e73e 549//____________________________________________________________________
3e778975 550Float_t* AliTRDseedV1::GetProbability(Bool_t force)
0906e73e 551{
3e778975 552 if(!force) return &fProb[0];
4d6aee34 553 if(!CookPID()) return NULL;
3e778975 554 return &fProb[0];
555}
556
557//____________________________________________________________
558Bool_t AliTRDseedV1::CookPID()
559{
0906e73e 560// Fill probability array for tracklet from the DB.
561//
562// Parameters
563//
564// Output
4d6aee34 565// returns pointer to the probability array and NULL if missing DB access
0906e73e 566//
2a3191bb 567// Retrieve PID probabilities for e+-, mu+-, K+-, pi+- and p+- from the DB according to tracklet information:
568// - estimated momentum at tracklet reference point
569// - dE/dx measurements
570// - tracklet length
571// - TRD layer
572// According to the steering settings specified in the reconstruction one of the following methods are used
573// - Neural Network [default] - option "nn"
574// - 2D Likelihood - option "!nn"
0906e73e 575
0906e73e 576 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
577 if (!calibration) {
578 AliError("No access to calibration data");
3e778975 579 return kFALSE;
0906e73e 580 }
581
4d6aee34 582 if (!fkReconstructor) {
3a039a31 583 AliError("Reconstructor not set.");
3e778975 584 return kFALSE;
4ba1d6ae 585 }
586
0906e73e 587 // Retrieve the CDB container class with the parametric detector response
4d6aee34 588 const AliTRDCalPID *pd = calibration->GetPIDObject(fkReconstructor->GetPIDMethod());
0906e73e 589 if (!pd) {
590 AliError("No access to AliTRDCalPID object");
3e778975 591 return kFALSE;
0906e73e 592 }
4d6aee34 593 //AliInfo(Form("Method[%d] : %s", fkReconstructor->GetRecoParam() ->GetPIDMethod(), pd->IsA()->GetName()));
10f75631 594
29b87567 595 // calculate tracklet length TO DO
0906e73e 596 Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
597 /// TMath::Sqrt((1.0 - fSnp[iPlane]*fSnp[iPlane]) / (1.0 + fTgl[iPlane]*fTgl[iPlane]));
598
599 //calculate dE/dx
4d6aee34 600 CookdEdx(fkReconstructor->GetNdEdxSlices());
0906e73e 601
602 // Sets the a priori probabilities
f83cd814 603 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++)
604 fProb[ispec] = pd->GetProbability(ispec, GetMomentum(), &fdEdx[0], length, GetPlane());
f301a656 605
3e778975 606 return kTRUE;
0906e73e 607}
608
609//____________________________________________________________________
e4f2f73d 610Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
611{
612 //
613 // Returns a quality measurement of the current seed
614 //
615
dd8059a8 616 Float_t zcorr = kZcorr ? GetTilt() * (fZfit[0] - fZref[0]) : 0.;
29b87567 617 return
3e778975 618 .5 * TMath::Abs(18.0 - GetN())
29b87567 619 + 10.* TMath::Abs(fYfit[1] - fYref[1])
620 + 5. * TMath::Abs(fYfit[0] - fYref[0] + zcorr)
dd8059a8 621 + 2. * TMath::Abs(fZfit[0] - fZref[0]) / GetPadLength();
e4f2f73d 622}
623
624//____________________________________________________________________
d937ad7a 625void AliTRDseedV1::GetCovAt(Double_t x, Double_t *cov) const
0906e73e 626{
d937ad7a 627// Computes covariance in the y-z plane at radial point x (in tracking coordinates)
628// and returns the results in the preallocated array cov[3] as :
629// cov[0] = Var(y)
630// cov[1] = Cov(yz)
631// cov[2] = Var(z)
632//
633// Details
634//
635// For the linear transformation
636// BEGIN_LATEX
637// Y = T_{x} X^{T}
638// END_LATEX
639// The error propagation has the general form
640// BEGIN_LATEX
641// C_{Y} = T_{x} C_{X} T_{x}^{T}
642// END_LATEX
643// We apply this formula 2 times. First to calculate the covariance of the tracklet
644// at point x we consider:
645// BEGIN_LATEX
646// T_{x} = (1 x); X=(y0 dy/dx); C_{X}=#(){#splitline{Var(y0) Cov(y0, dy/dx)}{Cov(y0, dy/dx) Var(dy/dx)}}
647// END_LATEX
648// and secondly to take into account the tilt angle
649// BEGIN_LATEX
650// T_{#alpha} = #(){#splitline{cos(#alpha) __ sin(#alpha)}{-sin(#alpha) __ cos(#alpha)}}; X=(y z); C_{X}=#(){#splitline{Var(y) 0}{0 Var(z)}}
651// END_LATEX
652//
653// using simple trigonometrics one can write for this last case
654// BEGIN_LATEX
655// 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})}}
656// END_LATEX
657// which can be aproximated for small alphas (2 deg) with
658// BEGIN_LATEX
659// 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}}}
660// END_LATEX
661//
662// before applying the tilt rotation we also apply systematic uncertainties to the tracklet
663// position which can be tunned from outside via the AliTRDrecoParam::SetSysCovMatrix(). They might
664// account for extra misalignment/miscalibration uncertainties.
665//
666// Author :
667// Alex Bercuci <A.Bercuci@gsi.de>
668// Date : Jan 8th 2009
669//
b1957d3c 670
671
d937ad7a 672 Double_t xr = fX0-x;
673 Double_t sy2 = fCov[0] +2.*xr*fCov[1] + xr*xr*fCov[2];
b72f4eaf 674 Double_t sz2 = fS2Z;
675 //GetPadLength()*GetPadLength()/12.;
0906e73e 676
d937ad7a 677 // insert systematic uncertainties
4d6aee34 678 if(fkReconstructor){
bb2db46c 679 Double_t sys[15]; memset(sys, 0, 15*sizeof(Double_t));
4d6aee34 680 fkReconstructor->GetRecoParam()->GetSysCovMatrix(sys);
bb2db46c 681 sy2 += sys[0];
682 sz2 += sys[1];
683 }
d937ad7a 684 // rotate covariance matrix
dd8059a8 685 Double_t t2 = GetTilt()*GetTilt();
d937ad7a 686 Double_t correction = 1./(1. + t2);
687 cov[0] = (sy2+t2*sz2)*correction;
dd8059a8 688 cov[1] = GetTilt()*(sz2 - sy2)*correction;
d937ad7a 689 cov[2] = (t2*sy2+sz2)*correction;
b72f4eaf 690
691 //printf("C(%6.1f %+6.3f %6.1f) [%s]\n", 1.e4*TMath::Sqrt(cov[0]), cov[1], 1.e4*TMath::Sqrt(cov[2]), IsRowCross()?" RC ":"-");
d937ad7a 692}
eb38ed55 693
bb2db46c 694//____________________________________________________________
4d6aee34 695Double_t AliTRDseedV1::GetCovSqrt(const Double_t * const c, Double_t *d)
bb2db46c 696{
697// Helper function to calculate the square root of the covariance matrix.
698// The input matrix is stored in the vector c and the result in the vector d.
41b7c7b6 699// Both arrays have to be initialized by the user with at least 3 elements. Return negative in case of failure.
bb2db46c 700//
ec3f0161 701// For calculating the square root of the symmetric matrix c
702// the following relation is used:
bb2db46c 703// BEGIN_LATEX
ec3f0161 704// C^{1/2} = VD^{1/2}V^{-1}
bb2db46c 705// END_LATEX
41b7c7b6 706// with V being the matrix with the n eigenvectors as columns.
ec3f0161 707// In case C is symmetric the followings are true:
708// - matrix D is diagonal with the diagonal given by the eigenvalues of C
41b7c7b6 709// - V = V^{-1}
bb2db46c 710//
711// Author A.Bercuci <A.Bercuci@gsi.de>
712// Date Mar 19 2009
713
4d6aee34 714 Double_t l[2], // eigenvalues
715 v[3]; // eigenvectors
bb2db46c 716 // the secular equation and its solution :
717 // (c[0]-L)(c[2]-L)-c[1]^2 = 0
718 // L^2 - L*Tr(c)+DET(c) = 0
719 // L12 = [Tr(c) +- sqrt(Tr(c)^2-4*DET(c))]/2
4d6aee34 720 Double_t tr = c[0]+c[2], // trace
721 det = c[0]*c[2]-c[1]*c[1]; // determinant
722 if(TMath::Abs(det)<1.e-20) return -1.;
723 Double_t dd = TMath::Sqrt(tr*tr - 4*det);
724 l[0] = .5*(tr + dd);
725 l[1] = .5*(tr - dd);
726 if(l[0]<0. || l[1]<0.) return -1.;
41b7c7b6 727
728 // the sym V matrix
729 // | v00 v10|
730 // | v10 v11|
4d6aee34 731 Double_t tmp = (l[0]-c[0])/c[1];
732 v[0] = TMath::Sqrt(1./(tmp*tmp+1));
733 v[1] = tmp*v[0];
734 v[2] = v[1]*c[1]/(l[1]-c[2]);
41b7c7b6 735 // the VD^{1/2}V is:
4d6aee34 736 l[0] = TMath::Sqrt(l[0]); l[1] = TMath::Sqrt(l[1]);
737 d[0] = v[0]*v[0]*l[0]+v[1]*v[1]*l[1];
738 d[1] = v[0]*v[1]*l[0]+v[1]*v[2]*l[1];
739 d[2] = v[1]*v[1]*l[0]+v[2]*v[2]*l[1];
bb2db46c 740
741 return 1.;
742}
743
744//____________________________________________________________
4d6aee34 745Double_t AliTRDseedV1::GetCovInv(const Double_t * const c, Double_t *d)
bb2db46c 746{
747// Helper function to calculate the inverse of the covariance matrix.
748// The input matrix is stored in the vector c and the result in the vector d.
749// Both arrays have to be initialized by the user with at least 3 elements
750// The return value is the determinant or 0 in case of singularity.
751//
752// Author A.Bercuci <A.Bercuci@gsi.de>
753// Date Mar 19 2009
754
4d6aee34 755 Double_t det = c[0]*c[2] - c[1]*c[1];
756 if(TMath::Abs(det)<1.e-20) return 0.;
757 Double_t invDet = 1./det;
758 d[0] = c[2]*invDet;
759 d[1] =-c[1]*invDet;
760 d[2] = c[0]*invDet;
761 return det;
bb2db46c 762}
0906e73e 763
d937ad7a 764//____________________________________________________________________
b72f4eaf 765UShort_t AliTRDseedV1::GetVolumeId() const
766{
767 Int_t ic=0;
768 while(ic<kNclusters && !fClusters[ic]) ic++;
769 return fClusters[ic] ? fClusters[ic]->GetVolumeId() : 0;
770}
771
f301a656 772//____________________________________________________________________
773TLinearFitter* AliTRDseedV1::GetFitterY()
774{
775 if(!fgFitterY) fgFitterY = new TLinearFitter(1, "pol1");
776 fgFitterY->ClearPoints();
777 return fgFitterY;
778}
779
780//____________________________________________________________________
781TLinearFitter* AliTRDseedV1::GetFitterZ()
782{
783 if(!fgFitterZ) fgFitterZ = new TLinearFitter(1, "pol1");
784 fgFitterZ->ClearPoints();
785 return fgFitterZ;
786}
b72f4eaf 787
788//____________________________________________________________________
e3cf3d02 789void AliTRDseedV1::Calibrate()
d937ad7a 790{
e3cf3d02 791// Retrieve calibration and position parameters from OCDB.
792// The following information are used
d937ad7a 793// - detector index
e3cf3d02 794// - column and row position of first attached cluster. If no clusters are attached
795// to the tracklet a random central chamber position (c=70, r=7) will be used.
796//
797// The following information is cached in the tracklet
798// t0 (trigger delay)
799// drift velocity
800// PRF width
801// omega*tau = tg(a_L)
802// diffusion coefficients (longitudinal and transversal)
d937ad7a 803//
804// Author :
805// Alex Bercuci <A.Bercuci@gsi.de>
806// Date : Jan 8th 2009
807//
eb38ed55 808
d937ad7a 809 AliCDBManager *cdb = AliCDBManager::Instance();
810 if(cdb->GetRun() < 0){
811 AliError("OCDB manager not properly initialized");
812 return;
813 }
0906e73e 814
e3cf3d02 815 AliTRDcalibDB *calib = AliTRDcalibDB::Instance();
816 AliTRDCalROC *vdROC = calib->GetVdriftROC(fDet),
817 *t0ROC = calib->GetT0ROC(fDet);;
818 const AliTRDCalDet *vdDet = calib->GetVdriftDet();
819 const AliTRDCalDet *t0Det = calib->GetT0Det();
d937ad7a 820
821 Int_t col = 70, row = 7;
822 AliTRDcluster **c = &fClusters[0];
3e778975 823 if(GetN()){
d937ad7a 824 Int_t ic = 0;
8d2bec9e 825 while (ic<kNclusters && !(*c)){ic++; c++;}
d937ad7a 826 if(*c){
827 col = (*c)->GetPadCol();
828 row = (*c)->GetPadRow();
829 }
830 }
3a039a31 831
e17f4785 832 fT0 = (t0Det->GetValue(fDet) + t0ROC->GetValue(col,row)) / AliTRDCommonParam::Instance()->GetSamplingFrequency();
e3cf3d02 833 fVD = vdDet->GetValue(fDet) * vdROC->GetValue(col, row);
834 fS2PRF = calib->GetPRFWidth(fDet, col, row); fS2PRF *= fS2PRF;
835 fExB = AliTRDCommonParam::Instance()->GetOmegaTau(fVD);
836 AliTRDCommonParam::Instance()->GetDiffCoeff(fDiffL,
837 fDiffT, fVD);
838 SetBit(kCalib, kTRUE);
0906e73e 839}
840
0906e73e 841//____________________________________________________________________
29b87567 842void AliTRDseedV1::SetOwner()
0906e73e 843{
29b87567 844 //AliInfo(Form("own [%s] fOwner[%s]", own?"YES":"NO", fOwner?"YES":"NO"));
845
846 if(TestBit(kOwner)) return;
8d2bec9e 847 for(int ic=0; ic<kNclusters; ic++){
29b87567 848 if(!fClusters[ic]) continue;
849 fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
850 }
851 SetBit(kOwner);
0906e73e 852}
853
eb2b4f91 854//____________________________________________________________
855void AliTRDseedV1::SetPadPlane(AliTRDpadPlane *p)
856{
857// Shortcut method to initialize pad geometry.
858 if(!p) return;
859 SetTilt(TMath::Tan(TMath::DegToRad()*p->GetTiltingAngle()));
860 SetPadLength(p->GetLengthIPad());
861 SetPadWidth(p->GetWidthIPad());
862}
863
864
e4f2f73d 865//____________________________________________________________________
4d6aee34 866Bool_t AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *const chamber, Bool_t tilt)
e4f2f73d 867{
1fd9389f 868//
869// Projective algorithm to attach clusters to seeding tracklets. The following steps are performed :
870// 1. Collapse x coordinate for the full detector plane
871// 2. truncated mean on y (r-phi) direction
872// 3. purge clusters
873// 4. truncated mean on z direction
874// 5. purge clusters
875//
876// Parameters
877// - chamber : pointer to tracking chamber container used to search the tracklet
878// - tilt : switch for tilt correction during road building [default true]
879// Output
880// - true : if tracklet found successfully. Failure can happend because of the following:
881// -
882// Detailed description
883//
884// We start up by defining the track direction in the xy plane and roads. The roads are calculated based
8a7ff53c 885// on tracking information (variance in the r-phi direction) and estimated variance of the standard
886// clusters (see AliTRDcluster::SetSigmaY2()) corrected for tilt (see GetCovAt()). From this the road is
887// BEGIN_LATEX
500851ab 888// r_{y} = 3*#sqrt{12*(#sigma^{2}_{Trk}(y) + #frac{#sigma^{2}_{cl}(y) + tg^{2}(#alpha_{L})#sigma^{2}_{cl}(z)}{1+tg^{2}(#alpha_{L})})}
8a7ff53c 889// r_{z} = 1.5*L_{pad}
890// END_LATEX
1fd9389f 891//
4b755889 892// Author : Alexandru Bercuci <A.Bercuci@gsi.de>
893// Debug : level >3
1fd9389f 894
4d6aee34 895 if(!fkReconstructor->GetRecoParam() ){
29b87567 896 AliError("Seed can not be used without a valid RecoParam.");
897 return kFALSE;
898 }
b1957d3c 899 // Initialize reco params for this tracklet
900 // 1. first time bin in the drift region
a2abcbc5 901 Int_t t0 = 14;
4d6aee34 902 Int_t kClmin = Int_t(fkReconstructor->GetRecoParam() ->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
29b87567 903
4d6aee34 904 Double_t sysCov[5]; fkReconstructor->GetRecoParam()->GetSysCovMatrix(sysCov);
8a7ff53c 905 Double_t s2yTrk= fRefCov[0],
906 s2yCl = 0.,
907 s2zCl = GetPadLength()*GetPadLength()/12.,
908 syRef = TMath::Sqrt(s2yTrk),
909 t2 = GetTilt()*GetTilt();
29b87567 910 //define roads
4d6aee34 911 Double_t kroady = 1., //fkReconstructor->GetRecoParam() ->GetRoad1y();
566bf887 912 kroadz = GetPadLength() * fkReconstructor->GetRecoParam()->GetRoadzMultiplicator() + 1.;
8a7ff53c 913 // define probing cluster (the perfect cluster) and default calibration
914 Short_t sig[] = {0, 0, 10, 30, 10, 0,0};
915 AliTRDcluster cp(fDet, 6, 75, 0, sig, 0);
4d6aee34 916 if(fkReconstructor->IsHLT())cp.SetRPhiMethod(AliTRDcluster::kCOG);
8a7ff53c 917 Calibrate();
918
ee8fb199 919 AliDebug(4, "");
920 AliDebug(4, Form("syKalman[%f] rY[%f] rZ[%f]", syRef, kroady, kroadz));
29b87567 921
922 // working variables
b1957d3c 923 const Int_t kNrows = 16;
4b755889 924 const Int_t kNcls = 3*kNclusters; // buffer size
925 AliTRDcluster *clst[kNrows][kNcls];
3044dfe5 926 Bool_t blst[kNrows][kNcls];
4b755889 927 Double_t cond[4], dx, dy, yt, zt, yres[kNrows][kNcls];
928 Int_t idxs[kNrows][kNcls], ncl[kNrows], ncls = 0;
b1957d3c 929 memset(ncl, 0, kNrows*sizeof(Int_t));
4b755889 930 memset(yres, 0, kNrows*kNcls*sizeof(Double_t));
3044dfe5 931 memset(blst, 0, kNrows*kNcls*sizeof(Bool_t)); //this is 8 times faster to memset than "memset(clst, 0, kNrows*kNcls*sizeof(AliTRDcluster*))"
b1957d3c 932
29b87567 933 // Do cluster projection
4d6aee34 934 AliTRDcluster *c = NULL;
935 AliTRDchamberTimeBin *layer = NULL;
b1957d3c 936 Bool_t kBUFFER = kFALSE;
4b755889 937 for (Int_t it = 0; it < kNtb; it++) {
b1957d3c 938 if(!(layer = chamber->GetTB(it))) continue;
29b87567 939 if(!Int_t(*layer)) continue;
8a7ff53c 940 // get track projection at layers position
b1957d3c 941 dx = fX0 - layer->GetX();
942 yt = fYref[0] - fYref[1] * dx;
943 zt = fZref[0] - fZref[1] * dx;
8a7ff53c 944 // get standard cluster error corrected for tilt
945 cp.SetLocalTimeBin(it);
946 cp.SetSigmaY2(0.02, fDiffT, fExB, dx, -1./*zt*/, fYref[1]);
d956a643 947 s2yCl = (cp.GetSigmaY2() + sysCov[0] + t2*s2zCl)/(1.+t2);
8a7ff53c 948 // get estimated road
949 kroady = 3.*TMath::Sqrt(12.*(s2yTrk + s2yCl));
950
ee8fb199 951 AliDebug(5, Form(" %2d x[%f] yt[%f] zt[%f]", it, dx, yt, zt));
952
953 AliDebug(5, Form(" syTrk[um]=%6.2f syCl[um]=%6.2f syClTlt[um]=%6.2f Ry[mm]=%f", 1.e4*TMath::Sqrt(s2yTrk), 1.e4*TMath::Sqrt(cp.GetSigmaY2()), 1.e4*TMath::Sqrt(s2yCl), 1.e1*kroady));
b1957d3c 954
8a7ff53c 955 // select clusters
b1957d3c 956 cond[0] = yt; cond[2] = kroady;
957 cond[1] = zt; cond[3] = kroadz;
958 Int_t n=0, idx[6];
959 layer->GetClusters(cond, idx, n, 6);
960 for(Int_t ic = n; ic--;){
961 c = (*layer)[idx[ic]];
962 dy = yt - c->GetY();
dd8059a8 963 dy += tilt ? GetTilt() * (c->GetZ() - zt) : 0.;
b1957d3c 964 // select clusters on a 3 sigmaKalman level
965/* if(tilt && TMath::Abs(dy) > 3.*syRef){
966 printf("too large !!!\n");
967 continue;
968 }*/
969 Int_t r = c->GetPadRow();
ee8fb199 970 AliDebug(5, Form(" -> dy[%f] yc[%f] r[%d]", TMath::Abs(dy), c->GetY(), r));
b1957d3c 971 clst[r][ncl[r]] = c;
3044dfe5 972 blst[r][ncl[r]] = kTRUE;
b1957d3c 973 idxs[r][ncl[r]] = idx[ic];
974 yres[r][ncl[r]] = dy;
975 ncl[r]++; ncls++;
976
4b755889 977 if(ncl[r] >= kNcls) {
978 AliWarning(Form("Cluster candidates reached buffer limit %d. Some may be lost.", kNcls));
b1957d3c 979 kBUFFER = kTRUE;
29b87567 980 break;
981 }
982 }
b1957d3c 983 if(kBUFFER) break;
29b87567 984 }
ee8fb199 985 AliDebug(4, Form("Found %d clusters. Processing ...", ncls));
986 if(ncls<kClmin){
987 AliDebug(2, Form("CLUSTERS FOUND %d LESS THAN THRESHOLD %d.", ncls, kClmin));
7c3eecb8 988 SetErrorMsg(kAttachClFound);
ee8fb199 989 return kFALSE;
990 }
991
b1957d3c 992 // analyze each row individualy
993 Double_t mean, syDis;
994 Int_t nrow[] = {0, 0, 0}, nr = 0, lr=-1;
995 for(Int_t ir=kNrows; ir--;){
996 if(!(ncl[ir])) continue;
997 if(lr>0 && lr-ir != 1){
ee8fb199 998 AliDebug(2, "Gap in rows attached");
29b87567 999 }
ee8fb199 1000 AliDebug(5, Form(" r[%d] n[%d]", ir, ncl[ir]));
b1957d3c 1001 // Evaluate truncated mean on the y direction
1002 if(ncl[ir] > 3) AliMathBase::EvaluateUni(ncl[ir], yres[ir], mean, syDis, Int_t(ncl[ir]*.8));
1003 else {
1004 mean = 0.; syDis = 0.;
8a7ff53c 1005 continue;
b1957d3c 1006 }
1007
a2fbb6ec 1008 if(fkReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) > 3 && fkReconstructor->IsDebugStreaming()){
1009 TTreeSRedirector &cstreamer = *fkReconstructor->GetDebugStream(AliTRDrecoParam::kTracker);
bb5265c9 1010 TVectorD vdy(ncl[ir], yres[ir]);
21a6e3ac 1011 UChar_t stat(0);
6ad5e6b2 1012 if(IsKink()) SETBIT(stat, 1);
1013 if(IsStandAlone()) SETBIT(stat, 2);
4b755889 1014 cstreamer << "AttachClusters"
21a6e3ac 1015 << "stat=" << stat
1016 << "det=" << fDet
1017 << "pt=" << fPt
1018 << "s2y=" << s2yTrk
bb5265c9 1019 << "dy=" << &vdy
21a6e3ac 1020 << "m=" << mean
1021 << "s=" << syDis
4b755889 1022 << "\n";
1023 }
1024
b1957d3c 1025 // TODO check mean and sigma agains cluster resolution !!
ee8fb199 1026 AliDebug(4, Form(" m[%f (%5.3fs)] s[%f]", mean, TMath::Abs(mean/syDis), syDis));
b1957d3c 1027 // select clusters on a 3 sigmaDistr level
1028 Bool_t kFOUND = kFALSE;
1029 for(Int_t ic = ncl[ir]; ic--;){
1030 if(yres[ir][ic] - mean > 3. * syDis){
3044dfe5 1031 blst[ir][ic] = kFALSE; continue;
b1957d3c 1032 }
1033 nrow[nr]++; kFOUND = kTRUE;
1034 }
1035 // exit loop
1036 if(kFOUND) nr++;
1037 lr = ir; if(nr>=3) break;
29b87567 1038 }
ee8fb199 1039 AliDebug(4, Form(" nr[%d = %d + %d + %d]", nr, nrow[0], nrow[1], nrow[2]));
b1957d3c 1040
1041 // classify cluster rows
1042 Int_t row = -1;
1043 switch(nr){
1044 case 1:
1045 row = lr;
1046 break;
1047 case 2:
1048 SetBit(kRowCross, kTRUE); // mark pad row crossing
1049 if(nrow[0] > nrow[1]){ row = lr+1; lr = -1;}
1050 else{
1051 row = lr; lr = 1;
1052 nrow[2] = nrow[1];
1053 nrow[1] = nrow[0];
1054 nrow[0] = nrow[2];
29b87567 1055 }
b1957d3c 1056 break;
1057 case 3:
1058 SetBit(kRowCross, kTRUE); // mark pad row crossing
1059 break;
29b87567 1060 }
ee8fb199 1061 AliDebug(4, Form(" Ncl[rowMax = %d] = %d", row, nrow[0]));
1062 if(row<0){
1063 AliDebug(2, Form("WRONG ROW %d.", row));
7c3eecb8 1064 SetErrorMsg(kAttachRow);
ee8fb199 1065 return kFALSE;
1066 }
b1957d3c 1067 // Select and store clusters
1068 // We should consider here :
1069 // 1. How far is the chamber boundary
1070 // 2. How big is the mean
3e778975 1071 Int_t n = 0;
b1957d3c 1072 for (Int_t ir = 0; ir < nr; ir++) {
1073 Int_t jr = row + ir*lr;
ee8fb199 1074 AliDebug(4, Form(" Ncl[%d] @ R[%d] attaching ...", ncl[jr], jr));
b1957d3c 1075 for (Int_t ic = 0; ic < ncl[jr]; ic++) {
3044dfe5 1076 if(!blst[jr][ic])continue;
1077 c = clst[jr][ic];
b1957d3c 1078 Int_t it = c->GetPadTime();
6ad5e6b2 1079 Int_t idx = it+kNtb*ir;
1080 if(fClusters[idx]){
1081 AliDebug(2, Form("Cluster position already allocated tb[%2d] r[%d]. Skip !", it, jr));
1082 SetErrorMsg(kAttachMultipleCl);
1083 continue;
1084 }
1085
b1957d3c 1086 // TODO proper indexing of clusters !!
6ad5e6b2 1087 fIndexes[idx] = chamber->GetTB(it)->GetGlobalIndex(idxs[jr][ic]);
1088 fClusters[idx] = c;
3e778975 1089 n++;
b1957d3c 1090 }
1091 }
6ad5e6b2 1092 SetN(n);
b1957d3c 1093
29b87567 1094 // number of minimum numbers of clusters expected for the tracklet
6ad5e6b2 1095 if (GetN() < kClmin){
1096 AliDebug(2, Form("NOT ENOUGH CLUSTERS ATTACHED TO THE TRACKLET %d [%d] FROM FOUND [%d].", GetN(), kClmin, n));
7c3eecb8 1097 SetErrorMsg(kAttachClAttach);
e4f2f73d 1098 return kFALSE;
1099 }
0906e73e 1100
e3cf3d02 1101 // Load calibration parameters for this tracklet
1102 Calibrate();
b1957d3c 1103
1104 // calculate dx for time bins in the drift region (calibration aware)
a2abcbc5 1105 Float_t x[2] = {0.,0.}; Int_t tb[2]={0,0};
1106 for (Int_t it = t0, irp=0; irp<2 && it < AliTRDtrackerV1::GetNTimeBins(); it++) {
b1957d3c 1107 if(!fClusters[it]) continue;
1108 x[irp] = fClusters[it]->GetX();
a2abcbc5 1109 tb[irp] = fClusters[it]->GetLocalTimeBin();
b1957d3c 1110 irp++;
e3cf3d02 1111 }
d86ed84c 1112 Int_t dtb = tb[1] - tb[0];
1113 fdX = dtb ? (x[0] - x[1]) / dtb : 0.15;
29b87567 1114 return kTRUE;
e4f2f73d 1115}
1116
03cef9b2 1117//____________________________________________________________
1118void AliTRDseedV1::Bootstrap(const AliTRDReconstructor *rec)
1119{
1120// Fill in all derived information. It has to be called after recovery from file or HLT.
1121// The primitive data are
1122// - list of clusters
1123// - detector (as the detector will be removed from clusters)
1124// - position of anode wire (fX0) - temporary
1125// - track reference position and direction
1126// - momentum of the track
1127// - time bin length [cm]
1128//
1129// A.Bercuci <A.Bercuci@gsi.de> Oct 30th 2008
1130//
4d6aee34 1131 fkReconstructor = rec;
03cef9b2 1132 AliTRDgeometry g;
1133 AliTRDpadPlane *pp = g.GetPadPlane(fDet);
dd8059a8 1134 fPad[0] = pp->GetLengthIPad();
1135 fPad[1] = pp->GetWidthIPad();
1136 fPad[3] = TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle());
e3cf3d02 1137 //fSnp = fYref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);
1138 //fTgl = fZref[1];
3e778975 1139 Int_t n = 0, nshare = 0, nused = 0;
03cef9b2 1140 AliTRDcluster **cit = &fClusters[0];
8d2bec9e 1141 for(Int_t ic = kNclusters; ic--; cit++){
03cef9b2 1142 if(!(*cit)) return;
3e778975 1143 n++;
1144 if((*cit)->IsShared()) nshare++;
1145 if((*cit)->IsUsed()) nused++;
03cef9b2 1146 }
3e778975 1147 SetN(n); SetNUsed(nused); SetNShared(nshare);
e3cf3d02 1148 Fit();
03cef9b2 1149 CookLabels();
1150 GetProbability();
1151}
1152
1153
e4f2f73d 1154//____________________________________________________________________
b72f4eaf 1155Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr)
e4f2f73d 1156{
16cca13f 1157//
1158// Linear fit of the clusters attached to the tracklet
1159//
1160// Parameters :
1161// - tilt : switch for tilt pad correction of cluster y position based on
1162// the z, dzdx info from outside [default false].
1163// - zcorr : switch for using z information to correct for anisochronity
1fd9389f 1164// and a finner error parameterization estimation [default false]
16cca13f 1165// Output :
1166// True if successful
1167//
1168// Detailed description
1169//
1170// Fit in the xy plane
1171//
1fd9389f 1172// The fit is performed to estimate the y position of the tracklet and the track
1173// angle in the bending plane. The clusters are represented in the chamber coordinate
1174// system (with respect to the anode wire - see AliTRDtrackerV1::FollowBackProlongation()
1175// on how this is set). The x and y position of the cluster and also their variances
1176// are known from clusterizer level (see AliTRDcluster::GetXloc(), AliTRDcluster::GetYloc(),
1177// AliTRDcluster::GetSX() and AliTRDcluster::GetSY()).
1178// If gaussian approximation is used to calculate y coordinate of the cluster the position
1179// is recalculated taking into account the track angle. The general formula to calculate the
1180// error of cluster position in the gaussian approximation taking into account diffusion and track
1181// inclination is given for TRD by:
1182// BEGIN_LATEX
1183// #sigma^{2}_{y} = #sigma^{2}_{PRF} + #frac{x#delta_{t}^{2}}{(1+tg(#alpha_{L}))^{2}} + #frac{x^{2}tg^{2}(#phi-#alpha_{L})tg^{2}(#alpha_{L})}{12}
1184// END_LATEX
1185//
1186// Since errors are calculated only in the y directions, radial errors (x direction) are mapped to y
1187// by projection i.e.
1188// BEGIN_LATEX
1189// #sigma_{x|y} = tg(#phi) #sigma_{x}
1190// END_LATEX
1191// and also by the lorentz angle correction
1192//
1193// Fit in the xz plane
1194//
1195// The "fit" is performed to estimate the radial position (x direction) where pad row cross happens.
1196// If no pad row crossing the z position is taken from geometry and radial position is taken from the xy
1197// fit (see below).
1198//
1199// There are two methods to estimate the radial position of the pad row cross:
1200// 1. leading cluster radial position : Here the lower part of the tracklet is considered and the last
1201// cluster registered (at radial x0) on this segment is chosen to mark the pad row crossing. The error
1202// of the z estimate is given by :
1203// BEGIN_LATEX
1204// #sigma_{z} = tg(#theta) #Delta x_{x_{0}}/12
1205// END_LATEX
1206// The systematic errors for this estimation are generated by the following sources:
1207// - no charge sharing between pad rows is considered (sharp cross)
1208// - missing cluster at row cross (noise peak-up, under-threshold signal etc.).
1209//
1210// 2. charge fit over the crossing point : Here the full energy deposit along the tracklet is considered
1211// to estimate the position of the crossing by a fit in the qx plane. The errors in the q directions are
1212// parameterized as s_q = q^2. The systematic errors for this estimation are generated by the following sources:
1213// - no general model for the qx dependence
1214// - physical fluctuations of the charge deposit
1215// - gain calibration dependence
1216//
1217// Estimation of the radial position of the tracklet
16cca13f 1218//
1fd9389f 1219// For pad row cross the radial position is taken from the xz fit (see above). Otherwise it is taken as the
1220// interpolation point of the tracklet i.e. the point where the error in y of the fit is minimum. The error
1221// in the y direction of the tracklet is (see AliTRDseedV1::GetCovAt()):
1222// BEGIN_LATEX
1223// #sigma_{y} = #sigma^{2}_{y_{0}} + 2xcov(y_{0}, dy/dx) + #sigma^{2}_{dy/dx}
1224// END_LATEX
1225// and thus the radial position is:
1226// BEGIN_LATEX
1227// x = - cov(y_{0}, dy/dx)/#sigma^{2}_{dy/dx}
1228// END_LATEX
1229//
1230// Estimation of tracklet position error
1231//
1232// The error in y direction is the error of the linear fit at the radial position of the tracklet while in the z
1233// direction is given by the cluster error or pad row cross error. In case of no pad row cross this is given by:
1234// BEGIN_LATEX
1235// #sigma_{y} = #sigma^{2}_{y_{0}} - 2cov^{2}(y_{0}, dy/dx)/#sigma^{2}_{dy/dx} + #sigma^{2}_{dy/dx}
1236// #sigma_{z} = Pad_{length}/12
1237// END_LATEX
1238// For pad row cross the full error is calculated at the radial position of the crossing (see above) and the error
1239// in z by the width of the crossing region - being a matter of parameterization.
1240// BEGIN_LATEX
1241// #sigma_{z} = tg(#theta) #Delta x_{x_{0}}/12
1242// END_LATEX
1243// In case of no tilt correction (default in the barrel tracking) the tilt is taken into account by the rotation of
1244// the covariance matrix. See AliTRDseedV1::GetCovAt() for details.
1245//
1246// Author
1247// A.Bercuci <A.Bercuci@gsi.de>
e4f2f73d 1248
b72f4eaf 1249 if(!IsCalibrated()) Calibrate();
e3cf3d02 1250
29b87567 1251 const Int_t kClmin = 8;
010d62b0 1252
2f7d6ac8 1253 // get track direction
1254 Double_t y0 = fYref[0];
1255 Double_t dydx = fYref[1];
1256 Double_t z0 = fZref[0];
1257 Double_t dzdx = fZref[1];
1258 Double_t yt, zt;
ae4e8b84 1259
5f1ae1e7 1260 AliTRDtrackerV1::AliTRDLeastSquare fitterY;
1261 AliTRDtrackerV1::AliTRDLeastSquare fitterZ;
f301a656 1262
29b87567 1263 // book cluster information
8d2bec9e 1264 Double_t qc[kNclusters], xc[kNclusters], yc[kNclusters], zc[kNclusters], sy[kNclusters];
e3cf3d02 1265
dd8059a8 1266 Int_t n = 0;
4d6aee34 1267 AliTRDcluster *c=NULL, **jc = &fClusters[0];
9eb2d46c 1268 for (Int_t ic=0; ic<kNtb; ic++, ++jc) {
29b87567 1269 xc[ic] = -1.;
1270 yc[ic] = 999.;
1271 zc[ic] = 999.;
1272 sy[ic] = 0.;
9eb2d46c 1273 if(!(c = (*jc))) continue;
29b87567 1274 if(!c->IsInChamber()) continue;
9462866a 1275
29b87567 1276 Float_t w = 1.;
1277 if(c->GetNPads()>4) w = .5;
1278 if(c->GetNPads()>5) w = .2;
010d62b0 1279
1fd9389f 1280 // cluster charge
dd8059a8 1281 qc[n] = TMath::Abs(c->GetQ());
1fd9389f 1282 // pad row of leading
1283
b72f4eaf 1284 // Radial cluster position
e3cf3d02 1285 //Int_t jc = TMath::Max(fN-3, 0);
1286 //xc[fN] = c->GetXloc(fT0, fVD, &qc[jc], &xc[jc]/*, z0 - c->GetX()*dzdx*/);
b72f4eaf 1287 xc[n] = fX0 - c->GetX();
1288
1fd9389f 1289 // extrapolated track to cluster position
dd8059a8 1290 yt = y0 - xc[n]*dydx;
dd8059a8 1291 zt = z0 - xc[n]*dzdx;
1fd9389f 1292
1293 // Recalculate cluster error based on tracking information
1294 c->SetSigmaY2(fS2PRF, fDiffT, fExB, xc[n], zcorr?zt:-1., dydx);
1295 sy[n] = TMath::Sqrt(c->GetSigmaY2());
1296
a2fbb6ec 1297 yc[n] = fkReconstructor->GetRecoParam()->UseGAUS() ?
1fd9389f 1298 c->GetYloc(y0, sy[n], GetPadWidth()): c->GetY();
1299 zc[n] = c->GetZ();
1300 //optional tilt correction
1301 if(tilt) yc[n] -= (GetTilt()*(zc[n] - zt));
1302
1303 fitterY.AddPoint(&xc[n], yc[n], TMath::Sqrt(sy[n]));
3044dfe5 1304 if(IsRowCross())fitterZ.AddPoint(&xc[n], qc[n], 1.);
dd8059a8 1305 n++;
29b87567 1306 }
3044dfe5 1307
47d5d320 1308 // to few clusters
dd8059a8 1309 if (n < kClmin) return kFALSE;
2f7d6ac8 1310
d937ad7a 1311 // fit XY
2f7d6ac8 1312 fitterY.Eval();
5f1ae1e7 1313 fYfit[0] = fitterY.GetFunctionParameter(0);
1314 fYfit[1] = -fitterY.GetFunctionParameter(1);
d937ad7a 1315 // store covariance
5f1ae1e7 1316 Double_t p[3];
1317 fitterY.GetCovarianceMatrix(p);
d937ad7a 1318 fCov[0] = p[0]; // variance of y0
5f1ae1e7 1319 fCov[1] = p[2]; // covariance of y0, dydx
1320 fCov[2] = p[1]; // variance of dydx
b1957d3c 1321 // the ref radial position is set at the minimum of
1322 // the y variance of the tracklet
b72f4eaf 1323 fX = -fCov[1]/fCov[2];
b1957d3c 1324
1325 // fit XZ
b72f4eaf 1326 if(IsRowCross()){
e355f67a 1327/* // THE LEADING CLUSTER METHOD
1fd9389f 1328 Float_t xMin = fX0;
b72f4eaf 1329 Int_t ic=n=kNclusters-1; jc = &fClusters[ic];
1fd9389f 1330 AliTRDcluster *c0 =0x0, **kc = &fClusters[kNtb-1];
1331 for(; ic>kNtb; ic--, --jc, --kc){
1332 if((c0 = (*kc)) && c0->IsInChamber() && (xMin>c0->GetX())) xMin = c0->GetX();
1333 if(!(c = (*jc))) continue;
1334 if(!c->IsInChamber()) continue;
1335 zc[kNclusters-1] = c->GetZ();
1336 fX = fX0 - c->GetX();
1337 }
1338 fZfit[0] = .5*(zc[0]+zc[kNclusters-1]); fZfit[1] = 0.;
1339 // Error parameterization
1340 fS2Z = fdX*fZref[1];
e355f67a 1341 fS2Z *= fS2Z; fS2Z *= 0.2887; // 1/sqrt(12)*/
1342
1fd9389f 1343 // THE FIT X-Q PLANE METHOD
e355f67a 1344 Int_t ic=n=kNclusters-1; jc = &fClusters[ic];
b72f4eaf 1345 for(; ic>kNtb; ic--, --jc){
1346 if(!(c = (*jc))) continue;
1347 if(!c->IsInChamber()) continue;
1348 qc[n] = TMath::Abs(c->GetQ());
1349 xc[n] = fX0 - c->GetX();
1350 zc[n] = c->GetZ();
1351 fitterZ.AddPoint(&xc[n], -qc[n], 1.);
1352 n--;
1353 }
1354 // fit XZ
1355 fitterZ.Eval();
5f1ae1e7 1356 if(fitterZ.GetFunctionParameter(1)!=0.){
1357 fX = -fitterZ.GetFunctionParameter(0)/fitterZ.GetFunctionParameter(1);
b72f4eaf 1358 fX=(fX<0.)?0.:fX;
1359 Float_t dl = .5*AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght();
1360 fX=(fX> dl)?dl:fX;
07bbc13c 1361 fX-=.055; // TODO to be understood
b72f4eaf 1362 }
1363
1364 fZfit[0] = .5*(zc[0]+zc[kNclusters-1]); fZfit[1] = 0.;
c850c351 1365 // temporary external error parameterization
1366 fS2Z = 0.05+0.4*TMath::Abs(fZref[1]); fS2Z *= fS2Z;
1367 // TODO correct formula
1368 //fS2Z = sigma_x*TMath::Abs(fZref[1]);
b1957d3c 1369 } else {
1370 fZfit[0] = zc[0]; fZfit[1] = 0.;
dd8059a8 1371 fS2Z = GetPadLength()*GetPadLength()/12.;
29b87567 1372 }
b72f4eaf 1373 fS2Y = fCov[0] +2.*fX*fCov[1] + fX*fX*fCov[2];
29b87567 1374 return kTRUE;
e4f2f73d 1375}
1376
e4f2f73d 1377
f29f13a6 1378/*
e3cf3d02 1379//_____________________________________________________________________________
1380void AliTRDseedV1::FitMI()
1381{
1382//
1383// Fit the seed.
1384// Marian Ivanov's version
1385//
1386// linear fit on the y direction with respect to the reference direction.
1387// The residuals for each x (x = xc - x0) are deduced from:
1388// dy = y - yt (1)
1389// the tilting correction is written :
1390// y = yc + h*(zc-zt) (2)
1391// yt = y0+dy/dx*x (3)
1392// zt = z0+dz/dx*x (4)
1393// from (1),(2),(3) and (4)
1394// dy = yc - y0 - (dy/dx + h*dz/dx)*x + h*(zc-z0)
1395// the last term introduces the correction on y direction due to tilting pads. There are 2 ways to account for this:
1396// 1. use tilting correction for calculating the y
1397// 2. neglect tilting correction here and account for it in the error parametrization of the tracklet.
1398 const Float_t kRatio = 0.8;
1399 const Int_t kClmin = 5;
1400 const Float_t kmaxtan = 2;
1401
1402 if (TMath::Abs(fYref[1]) > kmaxtan){
1403 //printf("Exit: Abs(fYref[1]) = %3.3f, kmaxtan = %3.3f\n", TMath::Abs(fYref[1]), kmaxtan);
1404 return; // Track inclined too much
1405 }
1406
1407 Float_t sigmaexp = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction
dd8059a8 1408 Float_t ycrosscor = GetPadLength() * GetTilt() * 0.5; // Y correction for crossing
e3cf3d02 1409 Int_t fNChange = 0;
1410
1411 Double_t sumw;
1412 Double_t sumwx;
1413 Double_t sumwx2;
1414 Double_t sumwy;
1415 Double_t sumwxy;
1416 Double_t sumwz;
1417 Double_t sumwxz;
1418
1419 // Buffering: Leave it constant fot Performance issues
1420 Int_t zints[kNtb]; // Histograming of the z coordinate
1421 // Get 1 and second max probable coodinates in z
1422 Int_t zouts[2*kNtb];
1423 Float_t allowedz[kNtb]; // Allowed z for given time bin
1424 Float_t yres[kNtb]; // Residuals from reference
dd8059a8 1425 //Float_t anglecor = GetTilt() * fZref[1]; // Correction to the angle
e3cf3d02 1426
1427 Float_t pos[3*kNtb]; memset(pos, 0, 3*kNtb*sizeof(Float_t));
1428 Float_t *fX = &pos[0], *fY = &pos[kNtb], *fZ = &pos[2*kNtb];
1429
1430 Int_t fN = 0; AliTRDcluster *c = 0x0;
1431 fN2 = 0;
1432 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1433 yres[i] = 10000.0;
1434 if (!(c = fClusters[i])) continue;
1435 if(!c->IsInChamber()) continue;
1436 // Residual y
dd8059a8 1437 //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] + GetTilt()*(fZ[i] - fZref[0]);
e3cf3d02 1438 fX[i] = fX0 - c->GetX();
1439 fY[i] = c->GetY();
1440 fZ[i] = c->GetZ();
dd8059a8 1441 yres[i] = fY[i] - GetTilt()*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
e3cf3d02 1442 zints[fN] = Int_t(fZ[i]);
1443 fN++;
1444 }
1445
1446 if (fN < kClmin){
1447 //printf("Exit fN < kClmin: fN = %d\n", fN);
1448 return;
1449 }
1450 Int_t nz = AliTRDtrackerV1::Freq(fN, zints, zouts, kFALSE);
1451 Float_t fZProb = zouts[0];
1452 if (nz <= 1) zouts[3] = 0;
1453 if (zouts[1] + zouts[3] < kClmin) {
1454 //printf("Exit zouts[1] = %d, zouts[3] = %d\n",zouts[1],zouts[3]);
1455 return;
1456 }
1457
1458 // Z distance bigger than pad - length
1459 if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) zouts[3] = 0;
1460
1461 Int_t breaktime = -1;
1462 Bool_t mbefore = kFALSE;
1463 Int_t cumul[kNtb][2];
1464 Int_t counts[2] = { 0, 0 };
1465
1466 if (zouts[3] >= 3) {
1467
1468 //
1469 // Find the break time allowing one chage on pad-rows
1470 // with maximal number of accepted clusters
1471 //
1472 fNChange = 1;
1473 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1474 cumul[i][0] = counts[0];
1475 cumul[i][1] = counts[1];
1476 if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++;
1477 if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++;
1478 }
1479 Int_t maxcount = 0;
1480 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1481 Int_t after = cumul[AliTRDtrackerV1::GetNTimeBins()][0] - cumul[i][0];
1482 Int_t before = cumul[i][1];
1483 if (after + before > maxcount) {
1484 maxcount = after + before;
1485 breaktime = i;
1486 mbefore = kFALSE;
1487 }
1488 after = cumul[AliTRDtrackerV1::GetNTimeBins()-1][1] - cumul[i][1];
1489 before = cumul[i][0];
1490 if (after + before > maxcount) {
1491 maxcount = after + before;
1492 breaktime = i;
1493 mbefore = kTRUE;
1494 }
1495 }
1496 breaktime -= 1;
1497 }
1498
1499 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1500 if (i > breaktime) allowedz[i] = mbefore ? zouts[2] : zouts[0];
1501 if (i <= breaktime) allowedz[i] = (!mbefore) ? zouts[2] : zouts[0];
1502 }
1503
1504 if (((allowedz[0] > allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] < 0)) ||
1505 ((allowedz[0] < allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] > 0))) {
1506 //
1507 // Tracklet z-direction not in correspondance with track z direction
1508 //
1509 fNChange = 0;
1510 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1511 allowedz[i] = zouts[0]; // Only longest taken
1512 }
1513 }
1514
1515 if (fNChange > 0) {
1516 //
1517 // Cross pad -row tracklet - take the step change into account
1518 //
1519 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1520 if (!fClusters[i]) continue;
1521 if(!fClusters[i]->IsInChamber()) continue;
1522 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
1523 // Residual y
dd8059a8 1524 //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] + GetTilt()*(fZ[i] - fZref[0]);
1525 yres[i] = fY[i] - GetTilt()*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
f29f13a6 1526// if (TMath::Abs(fZ[i] - fZProb) > 2) {
dd8059a8 1527// if (fZ[i] > fZProb) yres[i] += GetTilt() * GetPadLength();
1528// if (fZ[i] < fZProb) yres[i] -= GetTilt() * GetPadLength();
f29f13a6 1529 }
e3cf3d02 1530 }
1531 }
1532
1533 Double_t yres2[kNtb];
1534 Double_t mean;
1535 Double_t sigma;
1536 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1537 if (!fClusters[i]) continue;
1538 if(!fClusters[i]->IsInChamber()) continue;
1539 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
1540 yres2[fN2] = yres[i];
1541 fN2++;
1542 }
1543 if (fN2 < kClmin) {
1544 //printf("Exit fN2 < kClmin: fN2 = %d\n", fN2);
1545 fN2 = 0;
1546 return;
1547 }
1548 AliMathBase::EvaluateUni(fN2,yres2,mean,sigma, Int_t(fN2*kRatio-2.));
1549 if (sigma < sigmaexp * 0.8) {
1550 sigma = sigmaexp;
1551 }
1552 //Float_t fSigmaY = sigma;
1553
1554 // Reset sums
1555 sumw = 0;
1556 sumwx = 0;
1557 sumwx2 = 0;
1558 sumwy = 0;
1559 sumwxy = 0;
1560 sumwz = 0;
1561 sumwxz = 0;
1562
1563 fN2 = 0;
1564 Float_t fMeanz = 0;
1565 Float_t fMPads = 0;
1566 fUsable = 0;
1567 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1568 if (!fClusters[i]) continue;
1569 if (!fClusters[i]->IsInChamber()) continue;
1570 if (TMath::Abs(fZ[i] - allowedz[i]) > 2){fClusters[i] = 0x0; continue;}
1571 if (TMath::Abs(yres[i] - mean) > 4.0 * sigma){fClusters[i] = 0x0; continue;}
1572 SETBIT(fUsable,i);
1573 fN2++;
1574 fMPads += fClusters[i]->GetNPads();
1575 Float_t weight = 1.0;
1576 if (fClusters[i]->GetNPads() > 4) weight = 0.5;
1577 if (fClusters[i]->GetNPads() > 5) weight = 0.2;
1578
1579
1580 Double_t x = fX[i];
1581 //printf("x = %7.3f dy = %7.3f fit %7.3f\n", x, yres[i], fY[i]-yres[i]);
1582
1583 sumw += weight;
1584 sumwx += x * weight;
1585 sumwx2 += x*x * weight;
1586 sumwy += weight * yres[i];
1587 sumwxy += weight * (yres[i]) * x;
1588 sumwz += weight * fZ[i];
1589 sumwxz += weight * fZ[i] * x;
1590
1591 }
1592
1593 if (fN2 < kClmin){
1594 //printf("Exit fN2 < kClmin(2): fN2 = %d\n",fN2);
1595 fN2 = 0;
1596 return;
1597 }
1598 fMeanz = sumwz / sumw;
1599 Float_t correction = 0;
1600 if (fNChange > 0) {
1601 // Tracklet on boundary
1602 if (fMeanz < fZProb) correction = ycrosscor;
1603 if (fMeanz > fZProb) correction = -ycrosscor;
1604 }
1605
1606 Double_t det = sumw * sumwx2 - sumwx * sumwx;
1607 fYfit[0] = (sumwx2 * sumwy - sumwx * sumwxy) / det;
1608 fYfit[1] = (sumw * sumwxy - sumwx * sumwy) / det;
1609
1610 fS2Y = 0;
1611 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1612 if (!TESTBIT(fUsable,i)) continue;
1613 Float_t delta = yres[i] - fYfit[0] - fYfit[1] * fX[i];
1614 fS2Y += delta*delta;
1615 }
1616 fS2Y = TMath::Sqrt(fS2Y / Float_t(fN2-2));
1617 // TEMPORARY UNTIL covariance properly calculated
1618 fS2Y = TMath::Max(fS2Y, Float_t(.1));
1619
1620 fZfit[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det;
1621 fZfit[1] = (sumw * sumwxz - sumwx * sumwz) / det;
1622// fYfitR[0] += fYref[0] + correction;
1623// fYfitR[1] += fYref[1];
1624// fYfit[0] = fYfitR[0];
1625 fYfit[1] = -fYfit[1];
1626
1627 UpdateUsed();
f29f13a6 1628}*/
e3cf3d02 1629
e4f2f73d 1630//___________________________________________________________________
203967fc 1631void AliTRDseedV1::Print(Option_t *o) const
e4f2f73d 1632{
1633 //
1634 // Printing the seedstatus
1635 //
1636
b72f4eaf 1637 AliInfo(Form("Det[%3d] X0[%7.2f] Pad{L[%5.2f] W[%5.2f] Tilt[%+6.2f]}", fDet, fX0, GetPadLength(), GetPadWidth(), GetTilt()));
dd8059a8 1638 AliInfo(Form("N[%2d] Nused[%2d] Nshared[%2d] [%d]", GetN(), GetNUsed(), GetNShared(), fN));
b72f4eaf 1639 AliInfo(Form("FLAGS : RC[%c] Kink[%c] SA[%c]", IsRowCross()?'y':'n', IsKink()?'y':'n', IsStandAlone()?'y':'n'));
dd8059a8 1640
1641 Double_t cov[3], x=GetX();
1642 GetCovAt(x, cov);
1643 AliInfo(" | x[cm] | y[cm] | z[cm] | dydx | dzdx |");
1644 AliInfo(Form("Fit | %7.2f | %7.2f+-%7.2f | %7.2f+-%7.2f| %5.2f | ----- |", x, GetY(), TMath::Sqrt(cov[0]), GetZ(), TMath::Sqrt(cov[2]), fYfit[1]));
16cca13f 1645 AliInfo(Form("Ref | %7.2f | %7.2f+-%7.2f | %7.2f+-%7.2f| %5.2f | %5.2f |", x, fYref[0]-fX*fYref[1], TMath::Sqrt(fRefCov[0]), fZref[0]-fX*fYref[1], TMath::Sqrt(fRefCov[2]), fYref[1], fZref[1]))
ee8fb199 1646 AliInfo(Form("P / Pt [GeV/c] = %f / %f", GetMomentum(), fPt));
1647 AliInfo(Form("dEdx [a.u.] = %f / %f / %f / %f / %f/ %f / %f / %f", fdEdx[0], fdEdx[1], fdEdx[2], fdEdx[3], fdEdx[4], fdEdx[5], fdEdx[6], fdEdx[7]));
1648 AliInfo(Form("PID = %5.3f / %5.3f / %5.3f / %5.3f / %5.3f", fProb[0], fProb[1], fProb[2], fProb[3], fProb[4]));
203967fc 1649
1650 if(strcmp(o, "a")!=0) return;
1651
4dc4dc2e 1652 AliTRDcluster* const* jc = &fClusters[0];
8d2bec9e 1653 for(int ic=0; ic<kNclusters; ic++, jc++) {
4dc4dc2e 1654 if(!(*jc)) continue;
203967fc 1655 (*jc)->Print(o);
4dc4dc2e 1656 }
e4f2f73d 1657}
47d5d320 1658
203967fc 1659
1660//___________________________________________________________________
1661Bool_t AliTRDseedV1::IsEqual(const TObject *o) const
1662{
1663 // Checks if current instance of the class has the same essential members
1664 // as the given one
1665
1666 if(!o) return kFALSE;
1667 const AliTRDseedV1 *inTracklet = dynamic_cast<const AliTRDseedV1*>(o);
1668 if(!inTracklet) return kFALSE;
1669
1670 for (Int_t i = 0; i < 2; i++){
e3cf3d02 1671 if ( fYref[i] != inTracklet->fYref[i] ) return kFALSE;
1672 if ( fZref[i] != inTracklet->fZref[i] ) return kFALSE;
203967fc 1673 }
1674
e3cf3d02 1675 if ( fS2Y != inTracklet->fS2Y ) return kFALSE;
dd8059a8 1676 if ( GetTilt() != inTracklet->GetTilt() ) return kFALSE;
1677 if ( GetPadLength() != inTracklet->GetPadLength() ) return kFALSE;
203967fc 1678
8d2bec9e 1679 for (Int_t i = 0; i < kNclusters; i++){
e3cf3d02 1680// if ( fX[i] != inTracklet->GetX(i) ) return kFALSE;
1681// if ( fY[i] != inTracklet->GetY(i) ) return kFALSE;
1682// if ( fZ[i] != inTracklet->GetZ(i) ) return kFALSE;
1683 if ( fIndexes[i] != inTracklet->fIndexes[i] ) return kFALSE;
203967fc 1684 }
f29f13a6 1685// if ( fUsable != inTracklet->fUsable ) return kFALSE;
203967fc 1686
1687 for (Int_t i=0; i < 2; i++){
e3cf3d02 1688 if ( fYfit[i] != inTracklet->fYfit[i] ) return kFALSE;
1689 if ( fZfit[i] != inTracklet->fZfit[i] ) return kFALSE;
1690 if ( fLabels[i] != inTracklet->fLabels[i] ) return kFALSE;
203967fc 1691 }
1692
e3cf3d02 1693/* if ( fMeanz != inTracklet->GetMeanz() ) return kFALSE;
1694 if ( fZProb != inTracklet->GetZProb() ) return kFALSE;*/
3e778975 1695 if ( fN != inTracklet->fN ) return kFALSE;
1696 //if ( fNUsed != inTracklet->fNUsed ) return kFALSE;
e3cf3d02 1697 //if ( fFreq != inTracklet->GetFreq() ) return kFALSE;
1698 //if ( fNChange != inTracklet->GetNChange() ) return kFALSE;
203967fc 1699
e3cf3d02 1700 if ( fC != inTracklet->fC ) return kFALSE;
1701 //if ( fCC != inTracklet->GetCC() ) return kFALSE;
1702 if ( fChi2 != inTracklet->fChi2 ) return kFALSE;
203967fc 1703 // if ( fChi2Z != inTracklet->GetChi2Z() ) return kFALSE;
1704
e3cf3d02 1705 if ( fDet != inTracklet->fDet ) return kFALSE;
b25a5e9e 1706 if ( fPt != inTracklet->fPt ) return kFALSE;
e3cf3d02 1707 if ( fdX != inTracklet->fdX ) return kFALSE;
203967fc 1708
8d2bec9e 1709 for (Int_t iCluster = 0; iCluster < kNclusters; iCluster++){
203967fc 1710 AliTRDcluster *curCluster = fClusters[iCluster];
e3cf3d02 1711 AliTRDcluster *inCluster = inTracklet->fClusters[iCluster];
203967fc 1712 if (curCluster && inCluster){
1713 if (! curCluster->IsEqual(inCluster) ) {
1714 curCluster->Print();
1715 inCluster->Print();
1716 return kFALSE;
1717 }
1718 } else {
1719 // if one cluster exists, and corresponding
1720 // in other tracklet doesn't - return kFALSE
1721 if(curCluster || inCluster) return kFALSE;
1722 }
1723 }
1724 return kTRUE;
1725}
5d401b45 1726