- changed track storage to AliHLTGlobalBarrelTrack
[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);
1012 if(IsKink()) SETBIT(stat, 0);
1013 if(IsStandAlone()) SETBIT(stat, 1);
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();
1079 // TODO proper indexing of clusters !!
e3cf3d02 1080 fIndexes[it+kNtb*ir] = chamber->GetTB(it)->GetGlobalIndex(idxs[jr][ic]);
1081 fClusters[it+kNtb*ir] = c;
29b87567 1082
b1957d3c 1083 //printf("\tid[%2d] it[%d] idx[%d]\n", ic, it, fIndexes[it]);
1084
3e778975 1085 n++;
b1957d3c 1086 }
1087 }
1088
29b87567 1089 // number of minimum numbers of clusters expected for the tracklet
3e778975 1090 if (n < kClmin){
ee8fb199 1091 AliDebug(2, Form("NOT ENOUGH CLUSTERS TO FIT THE TRACKLET %d [%d].", n, kClmin));
7c3eecb8 1092 SetErrorMsg(kAttachClAttach);
e4f2f73d 1093 return kFALSE;
1094 }
3e778975 1095 SetN(n);
0906e73e 1096
e3cf3d02 1097 // Load calibration parameters for this tracklet
1098 Calibrate();
b1957d3c 1099
1100 // calculate dx for time bins in the drift region (calibration aware)
a2abcbc5 1101 Float_t x[2] = {0.,0.}; Int_t tb[2]={0,0};
1102 for (Int_t it = t0, irp=0; irp<2 && it < AliTRDtrackerV1::GetNTimeBins(); it++) {
b1957d3c 1103 if(!fClusters[it]) continue;
1104 x[irp] = fClusters[it]->GetX();
a2abcbc5 1105 tb[irp] = fClusters[it]->GetLocalTimeBin();
b1957d3c 1106 irp++;
e3cf3d02 1107 }
d86ed84c 1108 Int_t dtb = tb[1] - tb[0];
1109 fdX = dtb ? (x[0] - x[1]) / dtb : 0.15;
29b87567 1110 return kTRUE;
e4f2f73d 1111}
1112
03cef9b2 1113//____________________________________________________________
1114void AliTRDseedV1::Bootstrap(const AliTRDReconstructor *rec)
1115{
1116// Fill in all derived information. It has to be called after recovery from file or HLT.
1117// The primitive data are
1118// - list of clusters
1119// - detector (as the detector will be removed from clusters)
1120// - position of anode wire (fX0) - temporary
1121// - track reference position and direction
1122// - momentum of the track
1123// - time bin length [cm]
1124//
1125// A.Bercuci <A.Bercuci@gsi.de> Oct 30th 2008
1126//
4d6aee34 1127 fkReconstructor = rec;
03cef9b2 1128 AliTRDgeometry g;
1129 AliTRDpadPlane *pp = g.GetPadPlane(fDet);
dd8059a8 1130 fPad[0] = pp->GetLengthIPad();
1131 fPad[1] = pp->GetWidthIPad();
1132 fPad[3] = TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle());
e3cf3d02 1133 //fSnp = fYref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);
1134 //fTgl = fZref[1];
3e778975 1135 Int_t n = 0, nshare = 0, nused = 0;
03cef9b2 1136 AliTRDcluster **cit = &fClusters[0];
8d2bec9e 1137 for(Int_t ic = kNclusters; ic--; cit++){
03cef9b2 1138 if(!(*cit)) return;
3e778975 1139 n++;
1140 if((*cit)->IsShared()) nshare++;
1141 if((*cit)->IsUsed()) nused++;
03cef9b2 1142 }
3e778975 1143 SetN(n); SetNUsed(nused); SetNShared(nshare);
e3cf3d02 1144 Fit();
03cef9b2 1145 CookLabels();
1146 GetProbability();
1147}
1148
1149
e4f2f73d 1150//____________________________________________________________________
b72f4eaf 1151Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr)
e4f2f73d 1152{
16cca13f 1153//
1154// Linear fit of the clusters attached to the tracklet
1155//
1156// Parameters :
1157// - tilt : switch for tilt pad correction of cluster y position based on
1158// the z, dzdx info from outside [default false].
1159// - zcorr : switch for using z information to correct for anisochronity
1fd9389f 1160// and a finner error parameterization estimation [default false]
16cca13f 1161// Output :
1162// True if successful
1163//
1164// Detailed description
1165//
1166// Fit in the xy plane
1167//
1fd9389f 1168// The fit is performed to estimate the y position of the tracklet and the track
1169// angle in the bending plane. The clusters are represented in the chamber coordinate
1170// system (with respect to the anode wire - see AliTRDtrackerV1::FollowBackProlongation()
1171// on how this is set). The x and y position of the cluster and also their variances
1172// are known from clusterizer level (see AliTRDcluster::GetXloc(), AliTRDcluster::GetYloc(),
1173// AliTRDcluster::GetSX() and AliTRDcluster::GetSY()).
1174// If gaussian approximation is used to calculate y coordinate of the cluster the position
1175// is recalculated taking into account the track angle. The general formula to calculate the
1176// error of cluster position in the gaussian approximation taking into account diffusion and track
1177// inclination is given for TRD by:
1178// BEGIN_LATEX
1179// #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}
1180// END_LATEX
1181//
1182// Since errors are calculated only in the y directions, radial errors (x direction) are mapped to y
1183// by projection i.e.
1184// BEGIN_LATEX
1185// #sigma_{x|y} = tg(#phi) #sigma_{x}
1186// END_LATEX
1187// and also by the lorentz angle correction
1188//
1189// Fit in the xz plane
1190//
1191// The "fit" is performed to estimate the radial position (x direction) where pad row cross happens.
1192// If no pad row crossing the z position is taken from geometry and radial position is taken from the xy
1193// fit (see below).
1194//
1195// There are two methods to estimate the radial position of the pad row cross:
1196// 1. leading cluster radial position : Here the lower part of the tracklet is considered and the last
1197// cluster registered (at radial x0) on this segment is chosen to mark the pad row crossing. The error
1198// of the z estimate is given by :
1199// BEGIN_LATEX
1200// #sigma_{z} = tg(#theta) #Delta x_{x_{0}}/12
1201// END_LATEX
1202// The systematic errors for this estimation are generated by the following sources:
1203// - no charge sharing between pad rows is considered (sharp cross)
1204// - missing cluster at row cross (noise peak-up, under-threshold signal etc.).
1205//
1206// 2. charge fit over the crossing point : Here the full energy deposit along the tracklet is considered
1207// to estimate the position of the crossing by a fit in the qx plane. The errors in the q directions are
1208// parameterized as s_q = q^2. The systematic errors for this estimation are generated by the following sources:
1209// - no general model for the qx dependence
1210// - physical fluctuations of the charge deposit
1211// - gain calibration dependence
1212//
1213// Estimation of the radial position of the tracklet
16cca13f 1214//
1fd9389f 1215// For pad row cross the radial position is taken from the xz fit (see above). Otherwise it is taken as the
1216// interpolation point of the tracklet i.e. the point where the error in y of the fit is minimum. The error
1217// in the y direction of the tracklet is (see AliTRDseedV1::GetCovAt()):
1218// BEGIN_LATEX
1219// #sigma_{y} = #sigma^{2}_{y_{0}} + 2xcov(y_{0}, dy/dx) + #sigma^{2}_{dy/dx}
1220// END_LATEX
1221// and thus the radial position is:
1222// BEGIN_LATEX
1223// x = - cov(y_{0}, dy/dx)/#sigma^{2}_{dy/dx}
1224// END_LATEX
1225//
1226// Estimation of tracklet position error
1227//
1228// The error in y direction is the error of the linear fit at the radial position of the tracklet while in the z
1229// direction is given by the cluster error or pad row cross error. In case of no pad row cross this is given by:
1230// BEGIN_LATEX
1231// #sigma_{y} = #sigma^{2}_{y_{0}} - 2cov^{2}(y_{0}, dy/dx)/#sigma^{2}_{dy/dx} + #sigma^{2}_{dy/dx}
1232// #sigma_{z} = Pad_{length}/12
1233// END_LATEX
1234// For pad row cross the full error is calculated at the radial position of the crossing (see above) and the error
1235// in z by the width of the crossing region - being a matter of parameterization.
1236// BEGIN_LATEX
1237// #sigma_{z} = tg(#theta) #Delta x_{x_{0}}/12
1238// END_LATEX
1239// In case of no tilt correction (default in the barrel tracking) the tilt is taken into account by the rotation of
1240// the covariance matrix. See AliTRDseedV1::GetCovAt() for details.
1241//
1242// Author
1243// A.Bercuci <A.Bercuci@gsi.de>
e4f2f73d 1244
b72f4eaf 1245 if(!IsCalibrated()) Calibrate();
e3cf3d02 1246
29b87567 1247 const Int_t kClmin = 8;
010d62b0 1248
2f7d6ac8 1249 // get track direction
1250 Double_t y0 = fYref[0];
1251 Double_t dydx = fYref[1];
1252 Double_t z0 = fZref[0];
1253 Double_t dzdx = fZref[1];
1254 Double_t yt, zt;
ae4e8b84 1255
5f1ae1e7 1256 AliTRDtrackerV1::AliTRDLeastSquare fitterY;
1257 AliTRDtrackerV1::AliTRDLeastSquare fitterZ;
f301a656 1258
29b87567 1259 // book cluster information
8d2bec9e 1260 Double_t qc[kNclusters], xc[kNclusters], yc[kNclusters], zc[kNclusters], sy[kNclusters];
e3cf3d02 1261
dd8059a8 1262 Int_t n = 0;
4d6aee34 1263 AliTRDcluster *c=NULL, **jc = &fClusters[0];
9eb2d46c 1264 for (Int_t ic=0; ic<kNtb; ic++, ++jc) {
29b87567 1265 xc[ic] = -1.;
1266 yc[ic] = 999.;
1267 zc[ic] = 999.;
1268 sy[ic] = 0.;
9eb2d46c 1269 if(!(c = (*jc))) continue;
29b87567 1270 if(!c->IsInChamber()) continue;
9462866a 1271
29b87567 1272 Float_t w = 1.;
1273 if(c->GetNPads()>4) w = .5;
1274 if(c->GetNPads()>5) w = .2;
010d62b0 1275
1fd9389f 1276 // cluster charge
dd8059a8 1277 qc[n] = TMath::Abs(c->GetQ());
1fd9389f 1278 // pad row of leading
1279
b72f4eaf 1280 // Radial cluster position
e3cf3d02 1281 //Int_t jc = TMath::Max(fN-3, 0);
1282 //xc[fN] = c->GetXloc(fT0, fVD, &qc[jc], &xc[jc]/*, z0 - c->GetX()*dzdx*/);
b72f4eaf 1283 xc[n] = fX0 - c->GetX();
1284
1fd9389f 1285 // extrapolated track to cluster position
dd8059a8 1286 yt = y0 - xc[n]*dydx;
dd8059a8 1287 zt = z0 - xc[n]*dzdx;
1fd9389f 1288
1289 // Recalculate cluster error based on tracking information
1290 c->SetSigmaY2(fS2PRF, fDiffT, fExB, xc[n], zcorr?zt:-1., dydx);
1291 sy[n] = TMath::Sqrt(c->GetSigmaY2());
1292
a2fbb6ec 1293 yc[n] = fkReconstructor->GetRecoParam()->UseGAUS() ?
1fd9389f 1294 c->GetYloc(y0, sy[n], GetPadWidth()): c->GetY();
1295 zc[n] = c->GetZ();
1296 //optional tilt correction
1297 if(tilt) yc[n] -= (GetTilt()*(zc[n] - zt));
1298
1299 fitterY.AddPoint(&xc[n], yc[n], TMath::Sqrt(sy[n]));
3044dfe5 1300 if(IsRowCross())fitterZ.AddPoint(&xc[n], qc[n], 1.);
dd8059a8 1301 n++;
29b87567 1302 }
3044dfe5 1303
47d5d320 1304 // to few clusters
dd8059a8 1305 if (n < kClmin) return kFALSE;
2f7d6ac8 1306
d937ad7a 1307 // fit XY
2f7d6ac8 1308 fitterY.Eval();
5f1ae1e7 1309 fYfit[0] = fitterY.GetFunctionParameter(0);
1310 fYfit[1] = -fitterY.GetFunctionParameter(1);
d937ad7a 1311 // store covariance
5f1ae1e7 1312 Double_t p[3];
1313 fitterY.GetCovarianceMatrix(p);
d937ad7a 1314 fCov[0] = p[0]; // variance of y0
5f1ae1e7 1315 fCov[1] = p[2]; // covariance of y0, dydx
1316 fCov[2] = p[1]; // variance of dydx
b1957d3c 1317 // the ref radial position is set at the minimum of
1318 // the y variance of the tracklet
b72f4eaf 1319 fX = -fCov[1]/fCov[2];
b1957d3c 1320
1321 // fit XZ
b72f4eaf 1322 if(IsRowCross()){
e355f67a 1323/* // THE LEADING CLUSTER METHOD
1fd9389f 1324 Float_t xMin = fX0;
b72f4eaf 1325 Int_t ic=n=kNclusters-1; jc = &fClusters[ic];
1fd9389f 1326 AliTRDcluster *c0 =0x0, **kc = &fClusters[kNtb-1];
1327 for(; ic>kNtb; ic--, --jc, --kc){
1328 if((c0 = (*kc)) && c0->IsInChamber() && (xMin>c0->GetX())) xMin = c0->GetX();
1329 if(!(c = (*jc))) continue;
1330 if(!c->IsInChamber()) continue;
1331 zc[kNclusters-1] = c->GetZ();
1332 fX = fX0 - c->GetX();
1333 }
1334 fZfit[0] = .5*(zc[0]+zc[kNclusters-1]); fZfit[1] = 0.;
1335 // Error parameterization
1336 fS2Z = fdX*fZref[1];
e355f67a 1337 fS2Z *= fS2Z; fS2Z *= 0.2887; // 1/sqrt(12)*/
1338
1fd9389f 1339 // THE FIT X-Q PLANE METHOD
e355f67a 1340 Int_t ic=n=kNclusters-1; jc = &fClusters[ic];
b72f4eaf 1341 for(; ic>kNtb; ic--, --jc){
1342 if(!(c = (*jc))) continue;
1343 if(!c->IsInChamber()) continue;
1344 qc[n] = TMath::Abs(c->GetQ());
1345 xc[n] = fX0 - c->GetX();
1346 zc[n] = c->GetZ();
1347 fitterZ.AddPoint(&xc[n], -qc[n], 1.);
1348 n--;
1349 }
1350 // fit XZ
1351 fitterZ.Eval();
5f1ae1e7 1352 if(fitterZ.GetFunctionParameter(1)!=0.){
1353 fX = -fitterZ.GetFunctionParameter(0)/fitterZ.GetFunctionParameter(1);
b72f4eaf 1354 fX=(fX<0.)?0.:fX;
1355 Float_t dl = .5*AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght();
1356 fX=(fX> dl)?dl:fX;
07bbc13c 1357 fX-=.055; // TODO to be understood
b72f4eaf 1358 }
1359
1360 fZfit[0] = .5*(zc[0]+zc[kNclusters-1]); fZfit[1] = 0.;
c850c351 1361 // temporary external error parameterization
1362 fS2Z = 0.05+0.4*TMath::Abs(fZref[1]); fS2Z *= fS2Z;
1363 // TODO correct formula
1364 //fS2Z = sigma_x*TMath::Abs(fZref[1]);
b1957d3c 1365 } else {
1366 fZfit[0] = zc[0]; fZfit[1] = 0.;
dd8059a8 1367 fS2Z = GetPadLength()*GetPadLength()/12.;
29b87567 1368 }
b72f4eaf 1369 fS2Y = fCov[0] +2.*fX*fCov[1] + fX*fX*fCov[2];
29b87567 1370 return kTRUE;
e4f2f73d 1371}
1372
e4f2f73d 1373
f29f13a6 1374/*
e3cf3d02 1375//_____________________________________________________________________________
1376void AliTRDseedV1::FitMI()
1377{
1378//
1379// Fit the seed.
1380// Marian Ivanov's version
1381//
1382// linear fit on the y direction with respect to the reference direction.
1383// The residuals for each x (x = xc - x0) are deduced from:
1384// dy = y - yt (1)
1385// the tilting correction is written :
1386// y = yc + h*(zc-zt) (2)
1387// yt = y0+dy/dx*x (3)
1388// zt = z0+dz/dx*x (4)
1389// from (1),(2),(3) and (4)
1390// dy = yc - y0 - (dy/dx + h*dz/dx)*x + h*(zc-z0)
1391// the last term introduces the correction on y direction due to tilting pads. There are 2 ways to account for this:
1392// 1. use tilting correction for calculating the y
1393// 2. neglect tilting correction here and account for it in the error parametrization of the tracklet.
1394 const Float_t kRatio = 0.8;
1395 const Int_t kClmin = 5;
1396 const Float_t kmaxtan = 2;
1397
1398 if (TMath::Abs(fYref[1]) > kmaxtan){
1399 //printf("Exit: Abs(fYref[1]) = %3.3f, kmaxtan = %3.3f\n", TMath::Abs(fYref[1]), kmaxtan);
1400 return; // Track inclined too much
1401 }
1402
1403 Float_t sigmaexp = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction
dd8059a8 1404 Float_t ycrosscor = GetPadLength() * GetTilt() * 0.5; // Y correction for crossing
e3cf3d02 1405 Int_t fNChange = 0;
1406
1407 Double_t sumw;
1408 Double_t sumwx;
1409 Double_t sumwx2;
1410 Double_t sumwy;
1411 Double_t sumwxy;
1412 Double_t sumwz;
1413 Double_t sumwxz;
1414
1415 // Buffering: Leave it constant fot Performance issues
1416 Int_t zints[kNtb]; // Histograming of the z coordinate
1417 // Get 1 and second max probable coodinates in z
1418 Int_t zouts[2*kNtb];
1419 Float_t allowedz[kNtb]; // Allowed z for given time bin
1420 Float_t yres[kNtb]; // Residuals from reference
dd8059a8 1421 //Float_t anglecor = GetTilt() * fZref[1]; // Correction to the angle
e3cf3d02 1422
1423 Float_t pos[3*kNtb]; memset(pos, 0, 3*kNtb*sizeof(Float_t));
1424 Float_t *fX = &pos[0], *fY = &pos[kNtb], *fZ = &pos[2*kNtb];
1425
1426 Int_t fN = 0; AliTRDcluster *c = 0x0;
1427 fN2 = 0;
1428 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1429 yres[i] = 10000.0;
1430 if (!(c = fClusters[i])) continue;
1431 if(!c->IsInChamber()) continue;
1432 // Residual y
dd8059a8 1433 //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] + GetTilt()*(fZ[i] - fZref[0]);
e3cf3d02 1434 fX[i] = fX0 - c->GetX();
1435 fY[i] = c->GetY();
1436 fZ[i] = c->GetZ();
dd8059a8 1437 yres[i] = fY[i] - GetTilt()*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
e3cf3d02 1438 zints[fN] = Int_t(fZ[i]);
1439 fN++;
1440 }
1441
1442 if (fN < kClmin){
1443 //printf("Exit fN < kClmin: fN = %d\n", fN);
1444 return;
1445 }
1446 Int_t nz = AliTRDtrackerV1::Freq(fN, zints, zouts, kFALSE);
1447 Float_t fZProb = zouts[0];
1448 if (nz <= 1) zouts[3] = 0;
1449 if (zouts[1] + zouts[3] < kClmin) {
1450 //printf("Exit zouts[1] = %d, zouts[3] = %d\n",zouts[1],zouts[3]);
1451 return;
1452 }
1453
1454 // Z distance bigger than pad - length
1455 if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) zouts[3] = 0;
1456
1457 Int_t breaktime = -1;
1458 Bool_t mbefore = kFALSE;
1459 Int_t cumul[kNtb][2];
1460 Int_t counts[2] = { 0, 0 };
1461
1462 if (zouts[3] >= 3) {
1463
1464 //
1465 // Find the break time allowing one chage on pad-rows
1466 // with maximal number of accepted clusters
1467 //
1468 fNChange = 1;
1469 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1470 cumul[i][0] = counts[0];
1471 cumul[i][1] = counts[1];
1472 if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++;
1473 if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++;
1474 }
1475 Int_t maxcount = 0;
1476 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1477 Int_t after = cumul[AliTRDtrackerV1::GetNTimeBins()][0] - cumul[i][0];
1478 Int_t before = cumul[i][1];
1479 if (after + before > maxcount) {
1480 maxcount = after + before;
1481 breaktime = i;
1482 mbefore = kFALSE;
1483 }
1484 after = cumul[AliTRDtrackerV1::GetNTimeBins()-1][1] - cumul[i][1];
1485 before = cumul[i][0];
1486 if (after + before > maxcount) {
1487 maxcount = after + before;
1488 breaktime = i;
1489 mbefore = kTRUE;
1490 }
1491 }
1492 breaktime -= 1;
1493 }
1494
1495 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1496 if (i > breaktime) allowedz[i] = mbefore ? zouts[2] : zouts[0];
1497 if (i <= breaktime) allowedz[i] = (!mbefore) ? zouts[2] : zouts[0];
1498 }
1499
1500 if (((allowedz[0] > allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] < 0)) ||
1501 ((allowedz[0] < allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] > 0))) {
1502 //
1503 // Tracklet z-direction not in correspondance with track z direction
1504 //
1505 fNChange = 0;
1506 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1507 allowedz[i] = zouts[0]; // Only longest taken
1508 }
1509 }
1510
1511 if (fNChange > 0) {
1512 //
1513 // Cross pad -row tracklet - take the step change into account
1514 //
1515 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1516 if (!fClusters[i]) continue;
1517 if(!fClusters[i]->IsInChamber()) continue;
1518 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
1519 // Residual y
dd8059a8 1520 //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] + GetTilt()*(fZ[i] - fZref[0]);
1521 yres[i] = fY[i] - GetTilt()*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
f29f13a6 1522// if (TMath::Abs(fZ[i] - fZProb) > 2) {
dd8059a8 1523// if (fZ[i] > fZProb) yres[i] += GetTilt() * GetPadLength();
1524// if (fZ[i] < fZProb) yres[i] -= GetTilt() * GetPadLength();
f29f13a6 1525 }
e3cf3d02 1526 }
1527 }
1528
1529 Double_t yres2[kNtb];
1530 Double_t mean;
1531 Double_t sigma;
1532 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1533 if (!fClusters[i]) continue;
1534 if(!fClusters[i]->IsInChamber()) continue;
1535 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
1536 yres2[fN2] = yres[i];
1537 fN2++;
1538 }
1539 if (fN2 < kClmin) {
1540 //printf("Exit fN2 < kClmin: fN2 = %d\n", fN2);
1541 fN2 = 0;
1542 return;
1543 }
1544 AliMathBase::EvaluateUni(fN2,yres2,mean,sigma, Int_t(fN2*kRatio-2.));
1545 if (sigma < sigmaexp * 0.8) {
1546 sigma = sigmaexp;
1547 }
1548 //Float_t fSigmaY = sigma;
1549
1550 // Reset sums
1551 sumw = 0;
1552 sumwx = 0;
1553 sumwx2 = 0;
1554 sumwy = 0;
1555 sumwxy = 0;
1556 sumwz = 0;
1557 sumwxz = 0;
1558
1559 fN2 = 0;
1560 Float_t fMeanz = 0;
1561 Float_t fMPads = 0;
1562 fUsable = 0;
1563 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1564 if (!fClusters[i]) continue;
1565 if (!fClusters[i]->IsInChamber()) continue;
1566 if (TMath::Abs(fZ[i] - allowedz[i]) > 2){fClusters[i] = 0x0; continue;}
1567 if (TMath::Abs(yres[i] - mean) > 4.0 * sigma){fClusters[i] = 0x0; continue;}
1568 SETBIT(fUsable,i);
1569 fN2++;
1570 fMPads += fClusters[i]->GetNPads();
1571 Float_t weight = 1.0;
1572 if (fClusters[i]->GetNPads() > 4) weight = 0.5;
1573 if (fClusters[i]->GetNPads() > 5) weight = 0.2;
1574
1575
1576 Double_t x = fX[i];
1577 //printf("x = %7.3f dy = %7.3f fit %7.3f\n", x, yres[i], fY[i]-yres[i]);
1578
1579 sumw += weight;
1580 sumwx += x * weight;
1581 sumwx2 += x*x * weight;
1582 sumwy += weight * yres[i];
1583 sumwxy += weight * (yres[i]) * x;
1584 sumwz += weight * fZ[i];
1585 sumwxz += weight * fZ[i] * x;
1586
1587 }
1588
1589 if (fN2 < kClmin){
1590 //printf("Exit fN2 < kClmin(2): fN2 = %d\n",fN2);
1591 fN2 = 0;
1592 return;
1593 }
1594 fMeanz = sumwz / sumw;
1595 Float_t correction = 0;
1596 if (fNChange > 0) {
1597 // Tracklet on boundary
1598 if (fMeanz < fZProb) correction = ycrosscor;
1599 if (fMeanz > fZProb) correction = -ycrosscor;
1600 }
1601
1602 Double_t det = sumw * sumwx2 - sumwx * sumwx;
1603 fYfit[0] = (sumwx2 * sumwy - sumwx * sumwxy) / det;
1604 fYfit[1] = (sumw * sumwxy - sumwx * sumwy) / det;
1605
1606 fS2Y = 0;
1607 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1608 if (!TESTBIT(fUsable,i)) continue;
1609 Float_t delta = yres[i] - fYfit[0] - fYfit[1] * fX[i];
1610 fS2Y += delta*delta;
1611 }
1612 fS2Y = TMath::Sqrt(fS2Y / Float_t(fN2-2));
1613 // TEMPORARY UNTIL covariance properly calculated
1614 fS2Y = TMath::Max(fS2Y, Float_t(.1));
1615
1616 fZfit[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det;
1617 fZfit[1] = (sumw * sumwxz - sumwx * sumwz) / det;
1618// fYfitR[0] += fYref[0] + correction;
1619// fYfitR[1] += fYref[1];
1620// fYfit[0] = fYfitR[0];
1621 fYfit[1] = -fYfit[1];
1622
1623 UpdateUsed();
f29f13a6 1624}*/
e3cf3d02 1625
e4f2f73d 1626//___________________________________________________________________
203967fc 1627void AliTRDseedV1::Print(Option_t *o) const
e4f2f73d 1628{
1629 //
1630 // Printing the seedstatus
1631 //
1632
b72f4eaf 1633 AliInfo(Form("Det[%3d] X0[%7.2f] Pad{L[%5.2f] W[%5.2f] Tilt[%+6.2f]}", fDet, fX0, GetPadLength(), GetPadWidth(), GetTilt()));
dd8059a8 1634 AliInfo(Form("N[%2d] Nused[%2d] Nshared[%2d] [%d]", GetN(), GetNUsed(), GetNShared(), fN));
b72f4eaf 1635 AliInfo(Form("FLAGS : RC[%c] Kink[%c] SA[%c]", IsRowCross()?'y':'n', IsKink()?'y':'n', IsStandAlone()?'y':'n'));
dd8059a8 1636
1637 Double_t cov[3], x=GetX();
1638 GetCovAt(x, cov);
1639 AliInfo(" | x[cm] | y[cm] | z[cm] | dydx | dzdx |");
1640 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 1641 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 1642 AliInfo(Form("P / Pt [GeV/c] = %f / %f", GetMomentum(), fPt));
1643 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]));
1644 AliInfo(Form("PID = %5.3f / %5.3f / %5.3f / %5.3f / %5.3f", fProb[0], fProb[1], fProb[2], fProb[3], fProb[4]));
203967fc 1645
1646 if(strcmp(o, "a")!=0) return;
1647
4dc4dc2e 1648 AliTRDcluster* const* jc = &fClusters[0];
8d2bec9e 1649 for(int ic=0; ic<kNclusters; ic++, jc++) {
4dc4dc2e 1650 if(!(*jc)) continue;
203967fc 1651 (*jc)->Print(o);
4dc4dc2e 1652 }
e4f2f73d 1653}
47d5d320 1654
203967fc 1655
1656//___________________________________________________________________
1657Bool_t AliTRDseedV1::IsEqual(const TObject *o) const
1658{
1659 // Checks if current instance of the class has the same essential members
1660 // as the given one
1661
1662 if(!o) return kFALSE;
1663 const AliTRDseedV1 *inTracklet = dynamic_cast<const AliTRDseedV1*>(o);
1664 if(!inTracklet) return kFALSE;
1665
1666 for (Int_t i = 0; i < 2; i++){
e3cf3d02 1667 if ( fYref[i] != inTracklet->fYref[i] ) return kFALSE;
1668 if ( fZref[i] != inTracklet->fZref[i] ) return kFALSE;
203967fc 1669 }
1670
e3cf3d02 1671 if ( fS2Y != inTracklet->fS2Y ) return kFALSE;
dd8059a8 1672 if ( GetTilt() != inTracklet->GetTilt() ) return kFALSE;
1673 if ( GetPadLength() != inTracklet->GetPadLength() ) return kFALSE;
203967fc 1674
8d2bec9e 1675 for (Int_t i = 0; i < kNclusters; i++){
e3cf3d02 1676// if ( fX[i] != inTracklet->GetX(i) ) return kFALSE;
1677// if ( fY[i] != inTracklet->GetY(i) ) return kFALSE;
1678// if ( fZ[i] != inTracklet->GetZ(i) ) return kFALSE;
1679 if ( fIndexes[i] != inTracklet->fIndexes[i] ) return kFALSE;
203967fc 1680 }
f29f13a6 1681// if ( fUsable != inTracklet->fUsable ) return kFALSE;
203967fc 1682
1683 for (Int_t i=0; i < 2; i++){
e3cf3d02 1684 if ( fYfit[i] != inTracklet->fYfit[i] ) return kFALSE;
1685 if ( fZfit[i] != inTracklet->fZfit[i] ) return kFALSE;
1686 if ( fLabels[i] != inTracklet->fLabels[i] ) return kFALSE;
203967fc 1687 }
1688
e3cf3d02 1689/* if ( fMeanz != inTracklet->GetMeanz() ) return kFALSE;
1690 if ( fZProb != inTracklet->GetZProb() ) return kFALSE;*/
3e778975 1691 if ( fN != inTracklet->fN ) return kFALSE;
1692 //if ( fNUsed != inTracklet->fNUsed ) return kFALSE;
e3cf3d02 1693 //if ( fFreq != inTracklet->GetFreq() ) return kFALSE;
1694 //if ( fNChange != inTracklet->GetNChange() ) return kFALSE;
203967fc 1695
e3cf3d02 1696 if ( fC != inTracklet->fC ) return kFALSE;
1697 //if ( fCC != inTracklet->GetCC() ) return kFALSE;
1698 if ( fChi2 != inTracklet->fChi2 ) return kFALSE;
203967fc 1699 // if ( fChi2Z != inTracklet->GetChi2Z() ) return kFALSE;
1700
e3cf3d02 1701 if ( fDet != inTracklet->fDet ) return kFALSE;
b25a5e9e 1702 if ( fPt != inTracklet->fPt ) return kFALSE;
e3cf3d02 1703 if ( fdX != inTracklet->fdX ) return kFALSE;
203967fc 1704
8d2bec9e 1705 for (Int_t iCluster = 0; iCluster < kNclusters; iCluster++){
203967fc 1706 AliTRDcluster *curCluster = fClusters[iCluster];
e3cf3d02 1707 AliTRDcluster *inCluster = inTracklet->fClusters[iCluster];
203967fc 1708 if (curCluster && inCluster){
1709 if (! curCluster->IsEqual(inCluster) ) {
1710 curCluster->Print();
1711 inCluster->Print();
1712 return kFALSE;
1713 }
1714 } else {
1715 // if one cluster exists, and corresponding
1716 // in other tracklet doesn't - return kFALSE
1717 if(curCluster || inCluster) return kFALSE;
1718 }
1719 }
1720 return kTRUE;
1721}
5d401b45 1722