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