Remove obsolete classes
[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
f301a656 599 if(fReconstructor->IsHLT()){
600 // this can be done here, because in HLT we have another NN
601 // don't run HLT with the normal NN because here we assume that the NN was trained with only two output neurons!
602 memset(fProb,0,AliPID::kSPECIES*sizeof(fProb[0]));
603 fProb[AliPID::kElectron] = pd->GetProbability(AliPID::kElectron, GetMomentum(), &fdEdx[0], length, GetPlane());
604 fProb[AliPID::kPion] = 1 - fProb[AliPID::kElectron];
0906e73e 605 }
f301a656 606 else
607 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++)
608 fProb[ispec] = pd->GetProbability(ispec, GetMomentum(), &fdEdx[0], length, GetPlane());
609
3e778975 610 return kTRUE;
0906e73e 611}
612
613//____________________________________________________________________
e4f2f73d 614Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
615{
616 //
617 // Returns a quality measurement of the current seed
618 //
619
dd8059a8 620 Float_t zcorr = kZcorr ? GetTilt() * (fZfit[0] - fZref[0]) : 0.;
29b87567 621 return
3e778975 622 .5 * TMath::Abs(18.0 - GetN())
29b87567 623 + 10.* TMath::Abs(fYfit[1] - fYref[1])
624 + 5. * TMath::Abs(fYfit[0] - fYref[0] + zcorr)
dd8059a8 625 + 2. * TMath::Abs(fZfit[0] - fZref[0]) / GetPadLength();
e4f2f73d 626}
627
628//____________________________________________________________________
d937ad7a 629void AliTRDseedV1::GetCovAt(Double_t x, Double_t *cov) const
0906e73e 630{
d937ad7a 631// Computes covariance in the y-z plane at radial point x (in tracking coordinates)
632// and returns the results in the preallocated array cov[3] as :
633// cov[0] = Var(y)
634// cov[1] = Cov(yz)
635// cov[2] = Var(z)
636//
637// Details
638//
639// For the linear transformation
640// BEGIN_LATEX
641// Y = T_{x} X^{T}
642// END_LATEX
643// The error propagation has the general form
644// BEGIN_LATEX
645// C_{Y} = T_{x} C_{X} T_{x}^{T}
646// END_LATEX
647// We apply this formula 2 times. First to calculate the covariance of the tracklet
648// at point x we consider:
649// BEGIN_LATEX
650// T_{x} = (1 x); X=(y0 dy/dx); C_{X}=#(){#splitline{Var(y0) Cov(y0, dy/dx)}{Cov(y0, dy/dx) Var(dy/dx)}}
651// END_LATEX
652// and secondly to take into account the tilt angle
653// BEGIN_LATEX
654// T_{#alpha} = #(){#splitline{cos(#alpha) __ sin(#alpha)}{-sin(#alpha) __ cos(#alpha)}}; X=(y z); C_{X}=#(){#splitline{Var(y) 0}{0 Var(z)}}
655// END_LATEX
656//
657// using simple trigonometrics one can write for this last case
658// BEGIN_LATEX
659// 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})}}
660// END_LATEX
661// which can be aproximated for small alphas (2 deg) with
662// BEGIN_LATEX
663// 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}}}
664// END_LATEX
665//
666// before applying the tilt rotation we also apply systematic uncertainties to the tracklet
667// position which can be tunned from outside via the AliTRDrecoParam::SetSysCovMatrix(). They might
668// account for extra misalignment/miscalibration uncertainties.
669//
670// Author :
671// Alex Bercuci <A.Bercuci@gsi.de>
672// Date : Jan 8th 2009
673//
b1957d3c 674
675
d937ad7a 676 Double_t xr = fX0-x;
677 Double_t sy2 = fCov[0] +2.*xr*fCov[1] + xr*xr*fCov[2];
b72f4eaf 678 Double_t sz2 = fS2Z;
679 //GetPadLength()*GetPadLength()/12.;
0906e73e 680
d937ad7a 681 // insert systematic uncertainties
bb2db46c 682 if(fReconstructor){
683 Double_t sys[15]; memset(sys, 0, 15*sizeof(Double_t));
684 fReconstructor->GetRecoParam()->GetSysCovMatrix(sys);
685 sy2 += sys[0];
686 sz2 += sys[1];
687 }
d937ad7a 688 // rotate covariance matrix
dd8059a8 689 Double_t t2 = GetTilt()*GetTilt();
d937ad7a 690 Double_t correction = 1./(1. + t2);
691 cov[0] = (sy2+t2*sz2)*correction;
dd8059a8 692 cov[1] = GetTilt()*(sz2 - sy2)*correction;
d937ad7a 693 cov[2] = (t2*sy2+sz2)*correction;
b72f4eaf 694
695 //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 696}
eb38ed55 697
bb2db46c 698//____________________________________________________________
699Double_t AliTRDseedV1::GetCovSqrt(Double_t *c, Double_t *d)
700{
701// Helper function to calculate the square root of the covariance matrix.
702// The input matrix is stored in the vector c and the result in the vector d.
41b7c7b6 703// Both arrays have to be initialized by the user with at least 3 elements. Return negative in case of failure.
bb2db46c 704//
ec3f0161 705// For calculating the square root of the symmetric matrix c
706// the following relation is used:
bb2db46c 707// BEGIN_LATEX
ec3f0161 708// C^{1/2} = VD^{1/2}V^{-1}
bb2db46c 709// END_LATEX
41b7c7b6 710// with V being the matrix with the n eigenvectors as columns.
ec3f0161 711// In case C is symmetric the followings are true:
712// - matrix D is diagonal with the diagonal given by the eigenvalues of C
41b7c7b6 713// - V = V^{-1}
bb2db46c 714//
715// Author A.Bercuci <A.Bercuci@gsi.de>
716// Date Mar 19 2009
717
41b7c7b6 718 Double_t L[2], // eigenvalues
719 V[3]; // eigenvectors
bb2db46c 720 // the secular equation and its solution :
721 // (c[0]-L)(c[2]-L)-c[1]^2 = 0
722 // L^2 - L*Tr(c)+DET(c) = 0
723 // L12 = [Tr(c) +- sqrt(Tr(c)^2-4*DET(c))]/2
724 Double_t Tr = c[0]+c[2], // trace
725 DET = c[0]*c[2]-c[1]*c[1]; // determinant
41b7c7b6 726 if(TMath::Abs(DET)<1.e-20) return -1.;
bb2db46c 727 Double_t DD = TMath::Sqrt(Tr*Tr - 4*DET);
728 L[0] = .5*(Tr + DD);
729 L[1] = .5*(Tr - DD);
41b7c7b6 730 if(L[0]<0. || L[1]<0.) return -1.;
731
732 // the sym V matrix
733 // | v00 v10|
734 // | v10 v11|
735 Double_t tmp = (L[0]-c[0])/c[1];
736 V[0] = TMath::Sqrt(1./(tmp*tmp+1));
737 V[1] = tmp*V[0];
738 V[2] = V[1]*c[1]/(L[1]-c[2]);
739 // the VD^{1/2}V is:
740 L[0] = TMath::Sqrt(L[0]); L[1] = TMath::Sqrt(L[1]);
741 d[0] = V[0]*V[0]*L[0]+V[1]*V[1]*L[1];
742 d[1] = V[0]*V[1]*L[0]+V[1]*V[2]*L[1];
743 d[2] = V[1]*V[1]*L[0]+V[2]*V[2]*L[1];
bb2db46c 744
745 return 1.;
746}
747
748//____________________________________________________________
749Double_t AliTRDseedV1::GetCovInv(Double_t *c, Double_t *d)
750{
751// Helper function to calculate the inverse of the covariance matrix.
752// The input matrix is stored in the vector c and the result in the vector d.
753// Both arrays have to be initialized by the user with at least 3 elements
754// The return value is the determinant or 0 in case of singularity.
755//
756// Author A.Bercuci <A.Bercuci@gsi.de>
757// Date Mar 19 2009
758
759 Double_t Det = c[0]*c[2] - c[1]*c[1];
760 if(TMath::Abs(Det)<1.e-20) return 0.;
761 Double_t InvDet = 1./Det;
762 d[0] = c[2]*InvDet;
763 d[1] =-c[1]*InvDet;
764 d[2] = c[0]*InvDet;
765 return Det;
766}
0906e73e 767
d937ad7a 768//____________________________________________________________________
b72f4eaf 769UShort_t AliTRDseedV1::GetVolumeId() const
770{
771 Int_t ic=0;
772 while(ic<kNclusters && !fClusters[ic]) ic++;
773 return fClusters[ic] ? fClusters[ic]->GetVolumeId() : 0;
774}
775
f301a656 776//____________________________________________________________________
777TLinearFitter* AliTRDseedV1::GetFitterY()
778{
779 if(!fgFitterY) fgFitterY = new TLinearFitter(1, "pol1");
780 fgFitterY->ClearPoints();
781 return fgFitterY;
782}
783
784//____________________________________________________________________
785TLinearFitter* AliTRDseedV1::GetFitterZ()
786{
787 if(!fgFitterZ) fgFitterZ = new TLinearFitter(1, "pol1");
788 fgFitterZ->ClearPoints();
789 return fgFitterZ;
790}
b72f4eaf 791
792//____________________________________________________________________
e3cf3d02 793void AliTRDseedV1::Calibrate()
d937ad7a 794{
e3cf3d02 795// Retrieve calibration and position parameters from OCDB.
796// The following information are used
d937ad7a 797// - detector index
e3cf3d02 798// - column and row position of first attached cluster. If no clusters are attached
799// to the tracklet a random central chamber position (c=70, r=7) will be used.
800//
801// The following information is cached in the tracklet
802// t0 (trigger delay)
803// drift velocity
804// PRF width
805// omega*tau = tg(a_L)
806// diffusion coefficients (longitudinal and transversal)
d937ad7a 807//
808// Author :
809// Alex Bercuci <A.Bercuci@gsi.de>
810// Date : Jan 8th 2009
811//
eb38ed55 812
d937ad7a 813 AliCDBManager *cdb = AliCDBManager::Instance();
814 if(cdb->GetRun() < 0){
815 AliError("OCDB manager not properly initialized");
816 return;
817 }
0906e73e 818
e3cf3d02 819 AliTRDcalibDB *calib = AliTRDcalibDB::Instance();
820 AliTRDCalROC *vdROC = calib->GetVdriftROC(fDet),
821 *t0ROC = calib->GetT0ROC(fDet);;
822 const AliTRDCalDet *vdDet = calib->GetVdriftDet();
823 const AliTRDCalDet *t0Det = calib->GetT0Det();
d937ad7a 824
825 Int_t col = 70, row = 7;
826 AliTRDcluster **c = &fClusters[0];
3e778975 827 if(GetN()){
d937ad7a 828 Int_t ic = 0;
8d2bec9e 829 while (ic<kNclusters && !(*c)){ic++; c++;}
d937ad7a 830 if(*c){
831 col = (*c)->GetPadCol();
832 row = (*c)->GetPadRow();
833 }
834 }
3a039a31 835
e3cf3d02 836 fT0 = t0Det->GetValue(fDet) + t0ROC->GetValue(col,row);
837 fVD = vdDet->GetValue(fDet) * vdROC->GetValue(col, row);
838 fS2PRF = calib->GetPRFWidth(fDet, col, row); fS2PRF *= fS2PRF;
839 fExB = AliTRDCommonParam::Instance()->GetOmegaTau(fVD);
840 AliTRDCommonParam::Instance()->GetDiffCoeff(fDiffL,
841 fDiffT, fVD);
842 SetBit(kCalib, kTRUE);
0906e73e 843}
844
0906e73e 845//____________________________________________________________________
29b87567 846void AliTRDseedV1::SetOwner()
0906e73e 847{
29b87567 848 //AliInfo(Form("own [%s] fOwner[%s]", own?"YES":"NO", fOwner?"YES":"NO"));
849
850 if(TestBit(kOwner)) return;
8d2bec9e 851 for(int ic=0; ic<kNclusters; ic++){
29b87567 852 if(!fClusters[ic]) continue;
853 fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
854 }
855 SetBit(kOwner);
0906e73e 856}
857
eb2b4f91 858//____________________________________________________________
859void AliTRDseedV1::SetPadPlane(AliTRDpadPlane *p)
860{
861// Shortcut method to initialize pad geometry.
862 if(!p) return;
863 SetTilt(TMath::Tan(TMath::DegToRad()*p->GetTiltingAngle()));
864 SetPadLength(p->GetLengthIPad());
865 SetPadWidth(p->GetWidthIPad());
866}
867
868
e4f2f73d 869//____________________________________________________________________
b1957d3c 870Bool_t AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *chamber, Bool_t tilt)
e4f2f73d 871{
1fd9389f 872//
873// Projective algorithm to attach clusters to seeding tracklets. The following steps are performed :
874// 1. Collapse x coordinate for the full detector plane
875// 2. truncated mean on y (r-phi) direction
876// 3. purge clusters
877// 4. truncated mean on z direction
878// 5. purge clusters
879//
880// Parameters
881// - chamber : pointer to tracking chamber container used to search the tracklet
882// - tilt : switch for tilt correction during road building [default true]
883// Output
884// - true : if tracklet found successfully. Failure can happend because of the following:
885// -
886// Detailed description
887//
888// We start up by defining the track direction in the xy plane and roads. The roads are calculated based
8a7ff53c 889// on tracking information (variance in the r-phi direction) and estimated variance of the standard
890// clusters (see AliTRDcluster::SetSigmaY2()) corrected for tilt (see GetCovAt()). From this the road is
891// BEGIN_LATEX
500851ab 892// 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 893// r_{z} = 1.5*L_{pad}
894// END_LATEX
1fd9389f 895//
4b755889 896// Author : Alexandru Bercuci <A.Bercuci@gsi.de>
897// Debug : level >3
1fd9389f 898
b1957d3c 899 Bool_t kPRINT = kFALSE;
29b87567 900 if(!fReconstructor->GetRecoParam() ){
901 AliError("Seed can not be used without a valid RecoParam.");
902 return kFALSE;
903 }
b1957d3c 904 // Initialize reco params for this tracklet
905 // 1. first time bin in the drift region
a2abcbc5 906 Int_t t0 = 14;
b1957d3c 907 Int_t kClmin = Int_t(fReconstructor->GetRecoParam() ->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
29b87567 908
6e39bde4 909 Double_t sysCov[5]; fReconstructor->GetRecoParam()->GetSysCovMatrix(sysCov);
8a7ff53c 910 Double_t s2yTrk= fRefCov[0],
911 s2yCl = 0.,
912 s2zCl = GetPadLength()*GetPadLength()/12.,
913 syRef = TMath::Sqrt(s2yTrk),
914 t2 = GetTilt()*GetTilt();
29b87567 915 //define roads
8a7ff53c 916 Double_t kroady = 1., //fReconstructor->GetRecoParam() ->GetRoad1y();
917 kroadz = GetPadLength() * 1.5 + 1.;
918 // define probing cluster (the perfect cluster) and default calibration
919 Short_t sig[] = {0, 0, 10, 30, 10, 0,0};
920 AliTRDcluster cp(fDet, 6, 75, 0, sig, 0);
d07eec70 921 if(fReconstructor->IsHLT())cp.SetRPhiMethod(AliTRDcluster::kCOG);
8a7ff53c 922 Calibrate();
923
b1957d3c 924 if(kPRINT) printf("AttachClusters() sy[%f] road[%f]\n", syRef, kroady);
29b87567 925
926 // working variables
b1957d3c 927 const Int_t kNrows = 16;
4b755889 928 const Int_t kNcls = 3*kNclusters; // buffer size
929 AliTRDcluster *clst[kNrows][kNcls];
930 Double_t cond[4], dx, dy, yt, zt, yres[kNrows][kNcls];
931 Int_t idxs[kNrows][kNcls], ncl[kNrows], ncls = 0;
b1957d3c 932 memset(ncl, 0, kNrows*sizeof(Int_t));
4b755889 933 memset(yres, 0, kNrows*kNcls*sizeof(Double_t));
934 memset(clst, 0, kNrows*kNcls*sizeof(AliTRDcluster*));
b1957d3c 935
29b87567 936 // Do cluster projection
b1957d3c 937 AliTRDcluster *c = 0x0;
29b87567 938 AliTRDchamberTimeBin *layer = 0x0;
b1957d3c 939 Bool_t kBUFFER = kFALSE;
4b755889 940 for (Int_t it = 0; it < kNtb; it++) {
b1957d3c 941 if(!(layer = chamber->GetTB(it))) continue;
29b87567 942 if(!Int_t(*layer)) continue;
8a7ff53c 943 // get track projection at layers position
b1957d3c 944 dx = fX0 - layer->GetX();
945 yt = fYref[0] - fYref[1] * dx;
946 zt = fZref[0] - fZref[1] * dx;
8a7ff53c 947 // get standard cluster error corrected for tilt
948 cp.SetLocalTimeBin(it);
949 cp.SetSigmaY2(0.02, fDiffT, fExB, dx, -1./*zt*/, fYref[1]);
6e39bde4 950 s2yCl = (cp.GetSigmaY2() + sysCov[0] + t2*s2zCl)/(1.+t2);
8a7ff53c 951 // get estimated road
952 kroady = 3.*TMath::Sqrt(12.*(s2yTrk + s2yCl));
953
954 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 955
8a7ff53c 956 // select clusters
b1957d3c 957 cond[0] = yt; cond[2] = kroady;
958 cond[1] = zt; cond[3] = kroadz;
959 Int_t n=0, idx[6];
960 layer->GetClusters(cond, idx, n, 6);
961 for(Int_t ic = n; ic--;){
962 c = (*layer)[idx[ic]];
963 dy = yt - c->GetY();
dd8059a8 964 dy += tilt ? GetTilt() * (c->GetZ() - zt) : 0.;
b1957d3c 965 // select clusters on a 3 sigmaKalman level
966/* if(tilt && TMath::Abs(dy) > 3.*syRef){
967 printf("too large !!!\n");
968 continue;
969 }*/
970 Int_t r = c->GetPadRow();
971 if(kPRINT) printf("\t\t%d dy[%f] yc[%f] r[%d]\n", ic, TMath::Abs(dy), c->GetY(), r);
972 clst[r][ncl[r]] = c;
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 }
b1957d3c 985 if(kPRINT) printf("Found %d clusters\n", ncls);
986 if(ncls<kClmin) return kFALSE;
987
988 // analyze each row individualy
989 Double_t mean, syDis;
990 Int_t nrow[] = {0, 0, 0}, nr = 0, lr=-1;
991 for(Int_t ir=kNrows; ir--;){
992 if(!(ncl[ir])) continue;
993 if(lr>0 && lr-ir != 1){
994 if(kPRINT) printf("W - gap in rows attached !!\n");
29b87567 995 }
b1957d3c 996 if(kPRINT) printf("\tir[%d] lr[%d] n[%d]\n", ir, lr, ncl[ir]);
997 // Evaluate truncated mean on the y direction
998 if(ncl[ir] > 3) AliMathBase::EvaluateUni(ncl[ir], yres[ir], mean, syDis, Int_t(ncl[ir]*.8));
999 else {
1000 mean = 0.; syDis = 0.;
8a7ff53c 1001 continue;
b1957d3c 1002 }
1003
4b755889 1004 if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 3){
1005 TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker);
bb5265c9 1006 TVectorD vdy(ncl[ir], yres[ir]);
21a6e3ac 1007 UChar_t stat(0);
1008 if(IsKink()) SETBIT(stat, 0);
1009 if(IsStandAlone()) SETBIT(stat, 1);
4b755889 1010 cstreamer << "AttachClusters"
21a6e3ac 1011 << "stat=" << stat
1012 << "det=" << fDet
1013 << "pt=" << fPt
1014 << "s2y=" << s2yTrk
bb5265c9 1015 << "dy=" << &vdy
21a6e3ac 1016 << "m=" << mean
1017 << "s=" << syDis
4b755889 1018 << "\n";
1019 }
1020
b1957d3c 1021 // TODO check mean and sigma agains cluster resolution !!
8a7ff53c 1022 if(kPRINT) printf("\tr[%2d] m[%f %5.3fsigma] s[%f]\n", ir, mean, TMath::Abs(mean/syDis), syDis);
b1957d3c 1023 // select clusters on a 3 sigmaDistr level
1024 Bool_t kFOUND = kFALSE;
1025 for(Int_t ic = ncl[ir]; ic--;){
1026 if(yres[ir][ic] - mean > 3. * syDis){
1027 clst[ir][ic] = 0x0; continue;
1028 }
1029 nrow[nr]++; kFOUND = kTRUE;
1030 }
1031 // exit loop
1032 if(kFOUND) nr++;
1033 lr = ir; if(nr>=3) break;
29b87567 1034 }
b1957d3c 1035 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]);
1036
1037 // classify cluster rows
1038 Int_t row = -1;
1039 switch(nr){
1040 case 1:
1041 row = lr;
1042 break;
1043 case 2:
1044 SetBit(kRowCross, kTRUE); // mark pad row crossing
1045 if(nrow[0] > nrow[1]){ row = lr+1; lr = -1;}
1046 else{
1047 row = lr; lr = 1;
1048 nrow[2] = nrow[1];
1049 nrow[1] = nrow[0];
1050 nrow[0] = nrow[2];
29b87567 1051 }
b1957d3c 1052 break;
1053 case 3:
1054 SetBit(kRowCross, kTRUE); // mark pad row crossing
1055 break;
29b87567 1056 }
b1957d3c 1057 if(kPRINT) printf("\trow[%d] n[%d]\n\n", row, nrow[0]);
1058 if(row<0) return kFALSE;
29b87567 1059
b1957d3c 1060 // Select and store clusters
1061 // We should consider here :
1062 // 1. How far is the chamber boundary
1063 // 2. How big is the mean
3e778975 1064 Int_t n = 0;
b1957d3c 1065 for (Int_t ir = 0; ir < nr; ir++) {
1066 Int_t jr = row + ir*lr;
1067 if(kPRINT) printf("\tattach %d clusters for row %d\n", ncl[jr], jr);
1068 for (Int_t ic = 0; ic < ncl[jr]; ic++) {
1069 if(!(c = clst[jr][ic])) continue;
1070 Int_t it = c->GetPadTime();
1071 // TODO proper indexing of clusters !!
e3cf3d02 1072 fIndexes[it+kNtb*ir] = chamber->GetTB(it)->GetGlobalIndex(idxs[jr][ic]);
1073 fClusters[it+kNtb*ir] = c;
29b87567 1074
b1957d3c 1075 //printf("\tid[%2d] it[%d] idx[%d]\n", ic, it, fIndexes[it]);
1076
3e778975 1077 n++;
b1957d3c 1078 }
1079 }
1080
29b87567 1081 // number of minimum numbers of clusters expected for the tracklet
3e778975 1082 if (n < kClmin){
1083 //AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", n, kClmin));
e4f2f73d 1084 return kFALSE;
1085 }
3e778975 1086 SetN(n);
0906e73e 1087
e3cf3d02 1088 // Load calibration parameters for this tracklet
1089 Calibrate();
b1957d3c 1090
1091 // calculate dx for time bins in the drift region (calibration aware)
a2abcbc5 1092 Float_t x[2] = {0.,0.}; Int_t tb[2]={0,0};
1093 for (Int_t it = t0, irp=0; irp<2 && it < AliTRDtrackerV1::GetNTimeBins(); it++) {
b1957d3c 1094 if(!fClusters[it]) continue;
1095 x[irp] = fClusters[it]->GetX();
a2abcbc5 1096 tb[irp] = fClusters[it]->GetLocalTimeBin();
b1957d3c 1097 irp++;
e3cf3d02 1098 }
d86ed84c 1099 Int_t dtb = tb[1] - tb[0];
1100 fdX = dtb ? (x[0] - x[1]) / dtb : 0.15;
29b87567 1101 return kTRUE;
e4f2f73d 1102}
1103
03cef9b2 1104//____________________________________________________________
1105void AliTRDseedV1::Bootstrap(const AliTRDReconstructor *rec)
1106{
1107// Fill in all derived information. It has to be called after recovery from file or HLT.
1108// The primitive data are
1109// - list of clusters
1110// - detector (as the detector will be removed from clusters)
1111// - position of anode wire (fX0) - temporary
1112// - track reference position and direction
1113// - momentum of the track
1114// - time bin length [cm]
1115//
1116// A.Bercuci <A.Bercuci@gsi.de> Oct 30th 2008
1117//
1118 fReconstructor = rec;
1119 AliTRDgeometry g;
1120 AliTRDpadPlane *pp = g.GetPadPlane(fDet);
dd8059a8 1121 fPad[0] = pp->GetLengthIPad();
1122 fPad[1] = pp->GetWidthIPad();
1123 fPad[3] = TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle());
e3cf3d02 1124 //fSnp = fYref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);
1125 //fTgl = fZref[1];
3e778975 1126 Int_t n = 0, nshare = 0, nused = 0;
03cef9b2 1127 AliTRDcluster **cit = &fClusters[0];
8d2bec9e 1128 for(Int_t ic = kNclusters; ic--; cit++){
03cef9b2 1129 if(!(*cit)) return;
3e778975 1130 n++;
1131 if((*cit)->IsShared()) nshare++;
1132 if((*cit)->IsUsed()) nused++;
03cef9b2 1133 }
3e778975 1134 SetN(n); SetNUsed(nused); SetNShared(nshare);
e3cf3d02 1135 Fit();
03cef9b2 1136 CookLabels();
1137 GetProbability();
1138}
1139
1140
e4f2f73d 1141//____________________________________________________________________
b72f4eaf 1142Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr)
e4f2f73d 1143{
16cca13f 1144//
1145// Linear fit of the clusters attached to the tracklet
1146//
1147// Parameters :
1148// - tilt : switch for tilt pad correction of cluster y position based on
1149// the z, dzdx info from outside [default false].
1150// - zcorr : switch for using z information to correct for anisochronity
1fd9389f 1151// and a finner error parameterization estimation [default false]
16cca13f 1152// Output :
1153// True if successful
1154//
1155// Detailed description
1156//
1157// Fit in the xy plane
1158//
1fd9389f 1159// The fit is performed to estimate the y position of the tracklet and the track
1160// angle in the bending plane. The clusters are represented in the chamber coordinate
1161// system (with respect to the anode wire - see AliTRDtrackerV1::FollowBackProlongation()
1162// on how this is set). The x and y position of the cluster and also their variances
1163// are known from clusterizer level (see AliTRDcluster::GetXloc(), AliTRDcluster::GetYloc(),
1164// AliTRDcluster::GetSX() and AliTRDcluster::GetSY()).
1165// If gaussian approximation is used to calculate y coordinate of the cluster the position
1166// is recalculated taking into account the track angle. The general formula to calculate the
1167// error of cluster position in the gaussian approximation taking into account diffusion and track
1168// inclination is given for TRD by:
1169// BEGIN_LATEX
1170// #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}
1171// END_LATEX
1172//
1173// Since errors are calculated only in the y directions, radial errors (x direction) are mapped to y
1174// by projection i.e.
1175// BEGIN_LATEX
1176// #sigma_{x|y} = tg(#phi) #sigma_{x}
1177// END_LATEX
1178// and also by the lorentz angle correction
1179//
1180// Fit in the xz plane
1181//
1182// The "fit" is performed to estimate the radial position (x direction) where pad row cross happens.
1183// If no pad row crossing the z position is taken from geometry and radial position is taken from the xy
1184// fit (see below).
1185//
1186// There are two methods to estimate the radial position of the pad row cross:
1187// 1. leading cluster radial position : Here the lower part of the tracklet is considered and the last
1188// cluster registered (at radial x0) on this segment is chosen to mark the pad row crossing. The error
1189// of the z estimate is given by :
1190// BEGIN_LATEX
1191// #sigma_{z} = tg(#theta) #Delta x_{x_{0}}/12
1192// END_LATEX
1193// The systematic errors for this estimation are generated by the following sources:
1194// - no charge sharing between pad rows is considered (sharp cross)
1195// - missing cluster at row cross (noise peak-up, under-threshold signal etc.).
1196//
1197// 2. charge fit over the crossing point : Here the full energy deposit along the tracklet is considered
1198// to estimate the position of the crossing by a fit in the qx plane. The errors in the q directions are
1199// parameterized as s_q = q^2. The systematic errors for this estimation are generated by the following sources:
1200// - no general model for the qx dependence
1201// - physical fluctuations of the charge deposit
1202// - gain calibration dependence
1203//
1204// Estimation of the radial position of the tracklet
16cca13f 1205//
1fd9389f 1206// For pad row cross the radial position is taken from the xz fit (see above). Otherwise it is taken as the
1207// interpolation point of the tracklet i.e. the point where the error in y of the fit is minimum. The error
1208// in the y direction of the tracklet is (see AliTRDseedV1::GetCovAt()):
1209// BEGIN_LATEX
1210// #sigma_{y} = #sigma^{2}_{y_{0}} + 2xcov(y_{0}, dy/dx) + #sigma^{2}_{dy/dx}
1211// END_LATEX
1212// and thus the radial position is:
1213// BEGIN_LATEX
1214// x = - cov(y_{0}, dy/dx)/#sigma^{2}_{dy/dx}
1215// END_LATEX
1216//
1217// Estimation of tracklet position error
1218//
1219// The error in y direction is the error of the linear fit at the radial position of the tracklet while in the z
1220// direction is given by the cluster error or pad row cross error. In case of no pad row cross this is given by:
1221// BEGIN_LATEX
1222// #sigma_{y} = #sigma^{2}_{y_{0}} - 2cov^{2}(y_{0}, dy/dx)/#sigma^{2}_{dy/dx} + #sigma^{2}_{dy/dx}
1223// #sigma_{z} = Pad_{length}/12
1224// END_LATEX
1225// For pad row cross the full error is calculated at the radial position of the crossing (see above) and the error
1226// in z by the width of the crossing region - being a matter of parameterization.
1227// BEGIN_LATEX
1228// #sigma_{z} = tg(#theta) #Delta x_{x_{0}}/12
1229// END_LATEX
1230// In case of no tilt correction (default in the barrel tracking) the tilt is taken into account by the rotation of
1231// the covariance matrix. See AliTRDseedV1::GetCovAt() for details.
1232//
1233// Author
1234// A.Bercuci <A.Bercuci@gsi.de>
e4f2f73d 1235
b72f4eaf 1236 if(!IsCalibrated()) Calibrate();
e3cf3d02 1237
29b87567 1238 const Int_t kClmin = 8;
010d62b0 1239
2f7d6ac8 1240 // get track direction
1241 Double_t y0 = fYref[0];
1242 Double_t dydx = fYref[1];
1243 Double_t z0 = fZref[0];
1244 Double_t dzdx = fZref[1];
1245 Double_t yt, zt;
ae4e8b84 1246
b1957d3c 1247 //AliTRDtrackerV1::AliTRDLeastSquare fitterZ;
f301a656 1248 TLinearFitter& fitterY=*GetFitterY();
1249 TLinearFitter& fitterZ=*GetFitterZ();
1250
29b87567 1251 // book cluster information
8d2bec9e 1252 Double_t qc[kNclusters], xc[kNclusters], yc[kNclusters], zc[kNclusters], sy[kNclusters];
e3cf3d02 1253
dd8059a8 1254 Int_t n = 0;
9eb2d46c 1255 AliTRDcluster *c=0x0, **jc = &fClusters[0];
9eb2d46c 1256 for (Int_t ic=0; ic<kNtb; ic++, ++jc) {
29b87567 1257 xc[ic] = -1.;
1258 yc[ic] = 999.;
1259 zc[ic] = 999.;
1260 sy[ic] = 0.;
9eb2d46c 1261 if(!(c = (*jc))) continue;
29b87567 1262 if(!c->IsInChamber()) continue;
9462866a 1263
29b87567 1264 Float_t w = 1.;
1265 if(c->GetNPads()>4) w = .5;
1266 if(c->GetNPads()>5) w = .2;
010d62b0 1267
1fd9389f 1268 // cluster charge
dd8059a8 1269 qc[n] = TMath::Abs(c->GetQ());
1fd9389f 1270 // pad row of leading
1271
b72f4eaf 1272 // Radial cluster position
e3cf3d02 1273 //Int_t jc = TMath::Max(fN-3, 0);
1274 //xc[fN] = c->GetXloc(fT0, fVD, &qc[jc], &xc[jc]/*, z0 - c->GetX()*dzdx*/);
b72f4eaf 1275 xc[n] = fX0 - c->GetX();
1276
1fd9389f 1277 // extrapolated track to cluster position
dd8059a8 1278 yt = y0 - xc[n]*dydx;
dd8059a8 1279 zt = z0 - xc[n]*dzdx;
1fd9389f 1280
1281 // Recalculate cluster error based on tracking information
1282 c->SetSigmaY2(fS2PRF, fDiffT, fExB, xc[n], zcorr?zt:-1., dydx);
1283 sy[n] = TMath::Sqrt(c->GetSigmaY2());
1284
1285 yc[n] = fReconstructor->UseGAUS() ?
1286 c->GetYloc(y0, sy[n], GetPadWidth()): c->GetY();
1287 zc[n] = c->GetZ();
1288 //optional tilt correction
1289 if(tilt) yc[n] -= (GetTilt()*(zc[n] - zt));
1290
1291 fitterY.AddPoint(&xc[n], yc[n], TMath::Sqrt(sy[n]));
b72f4eaf 1292 fitterZ.AddPoint(&xc[n], qc[n], 1.);
dd8059a8 1293 n++;
29b87567 1294 }
47d5d320 1295 // to few clusters
dd8059a8 1296 if (n < kClmin) return kFALSE;
2f7d6ac8 1297
d937ad7a 1298 // fit XY
2f7d6ac8 1299 fitterY.Eval();
d937ad7a 1300 fYfit[0] = fitterY.GetParameter(0);
1301 fYfit[1] = -fitterY.GetParameter(1);
1302 // store covariance
1303 Double_t *p = fitterY.GetCovarianceMatrix();
1304 fCov[0] = p[0]; // variance of y0
1305 fCov[1] = p[1]; // covariance of y0, dydx
1306 fCov[2] = p[3]; // variance of dydx
b1957d3c 1307 // the ref radial position is set at the minimum of
1308 // the y variance of the tracklet
b72f4eaf 1309 fX = -fCov[1]/fCov[2];
b1957d3c 1310
1311 // fit XZ
b72f4eaf 1312 if(IsRowCross()){
e355f67a 1313/* // THE LEADING CLUSTER METHOD
1fd9389f 1314 Float_t xMin = fX0;
b72f4eaf 1315 Int_t ic=n=kNclusters-1; jc = &fClusters[ic];
1fd9389f 1316 AliTRDcluster *c0 =0x0, **kc = &fClusters[kNtb-1];
1317 for(; ic>kNtb; ic--, --jc, --kc){
1318 if((c0 = (*kc)) && c0->IsInChamber() && (xMin>c0->GetX())) xMin = c0->GetX();
1319 if(!(c = (*jc))) continue;
1320 if(!c->IsInChamber()) continue;
1321 zc[kNclusters-1] = c->GetZ();
1322 fX = fX0 - c->GetX();
1323 }
1324 fZfit[0] = .5*(zc[0]+zc[kNclusters-1]); fZfit[1] = 0.;
1325 // Error parameterization
1326 fS2Z = fdX*fZref[1];
e355f67a 1327 fS2Z *= fS2Z; fS2Z *= 0.2887; // 1/sqrt(12)*/
1328
1fd9389f 1329 // THE FIT X-Q PLANE METHOD
e355f67a 1330 Int_t ic=n=kNclusters-1; jc = &fClusters[ic];
b72f4eaf 1331 for(; ic>kNtb; ic--, --jc){
1332 if(!(c = (*jc))) continue;
1333 if(!c->IsInChamber()) continue;
1334 qc[n] = TMath::Abs(c->GetQ());
1335 xc[n] = fX0 - c->GetX();
1336 zc[n] = c->GetZ();
1337 fitterZ.AddPoint(&xc[n], -qc[n], 1.);
1338 n--;
1339 }
1340 // fit XZ
1341 fitterZ.Eval();
1342 if(fitterZ.GetParameter(1)!=0.){
1343 fX = -fitterZ.GetParameter(0)/fitterZ.GetParameter(1);
1344 fX=(fX<0.)?0.:fX;
1345 Float_t dl = .5*AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght();
1346 fX=(fX> dl)?dl:fX;
07bbc13c 1347 fX-=.055; // TODO to be understood
b72f4eaf 1348 }
1349
1350 fZfit[0] = .5*(zc[0]+zc[kNclusters-1]); fZfit[1] = 0.;
c850c351 1351 // temporary external error parameterization
1352 fS2Z = 0.05+0.4*TMath::Abs(fZref[1]); fS2Z *= fS2Z;
1353 // TODO correct formula
1354 //fS2Z = sigma_x*TMath::Abs(fZref[1]);
b1957d3c 1355 } else {
1356 fZfit[0] = zc[0]; fZfit[1] = 0.;
dd8059a8 1357 fS2Z = GetPadLength()*GetPadLength()/12.;
29b87567 1358 }
b72f4eaf 1359 fS2Y = fCov[0] +2.*fX*fCov[1] + fX*fX*fCov[2];
29b87567 1360 return kTRUE;
e4f2f73d 1361}
1362
e4f2f73d 1363
f29f13a6 1364/*
e3cf3d02 1365//_____________________________________________________________________________
1366void AliTRDseedV1::FitMI()
1367{
1368//
1369// Fit the seed.
1370// Marian Ivanov's version
1371//
1372// linear fit on the y direction with respect to the reference direction.
1373// The residuals for each x (x = xc - x0) are deduced from:
1374// dy = y - yt (1)
1375// the tilting correction is written :
1376// y = yc + h*(zc-zt) (2)
1377// yt = y0+dy/dx*x (3)
1378// zt = z0+dz/dx*x (4)
1379// from (1),(2),(3) and (4)
1380// dy = yc - y0 - (dy/dx + h*dz/dx)*x + h*(zc-z0)
1381// the last term introduces the correction on y direction due to tilting pads. There are 2 ways to account for this:
1382// 1. use tilting correction for calculating the y
1383// 2. neglect tilting correction here and account for it in the error parametrization of the tracklet.
1384 const Float_t kRatio = 0.8;
1385 const Int_t kClmin = 5;
1386 const Float_t kmaxtan = 2;
1387
1388 if (TMath::Abs(fYref[1]) > kmaxtan){
1389 //printf("Exit: Abs(fYref[1]) = %3.3f, kmaxtan = %3.3f\n", TMath::Abs(fYref[1]), kmaxtan);
1390 return; // Track inclined too much
1391 }
1392
1393 Float_t sigmaexp = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction
dd8059a8 1394 Float_t ycrosscor = GetPadLength() * GetTilt() * 0.5; // Y correction for crossing
e3cf3d02 1395 Int_t fNChange = 0;
1396
1397 Double_t sumw;
1398 Double_t sumwx;
1399 Double_t sumwx2;
1400 Double_t sumwy;
1401 Double_t sumwxy;
1402 Double_t sumwz;
1403 Double_t sumwxz;
1404
1405 // Buffering: Leave it constant fot Performance issues
1406 Int_t zints[kNtb]; // Histograming of the z coordinate
1407 // Get 1 and second max probable coodinates in z
1408 Int_t zouts[2*kNtb];
1409 Float_t allowedz[kNtb]; // Allowed z for given time bin
1410 Float_t yres[kNtb]; // Residuals from reference
dd8059a8 1411 //Float_t anglecor = GetTilt() * fZref[1]; // Correction to the angle
e3cf3d02 1412
1413 Float_t pos[3*kNtb]; memset(pos, 0, 3*kNtb*sizeof(Float_t));
1414 Float_t *fX = &pos[0], *fY = &pos[kNtb], *fZ = &pos[2*kNtb];
1415
1416 Int_t fN = 0; AliTRDcluster *c = 0x0;
1417 fN2 = 0;
1418 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1419 yres[i] = 10000.0;
1420 if (!(c = fClusters[i])) continue;
1421 if(!c->IsInChamber()) continue;
1422 // Residual y
dd8059a8 1423 //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] + GetTilt()*(fZ[i] - fZref[0]);
e3cf3d02 1424 fX[i] = fX0 - c->GetX();
1425 fY[i] = c->GetY();
1426 fZ[i] = c->GetZ();
dd8059a8 1427 yres[i] = fY[i] - GetTilt()*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
e3cf3d02 1428 zints[fN] = Int_t(fZ[i]);
1429 fN++;
1430 }
1431
1432 if (fN < kClmin){
1433 //printf("Exit fN < kClmin: fN = %d\n", fN);
1434 return;
1435 }
1436 Int_t nz = AliTRDtrackerV1::Freq(fN, zints, zouts, kFALSE);
1437 Float_t fZProb = zouts[0];
1438 if (nz <= 1) zouts[3] = 0;
1439 if (zouts[1] + zouts[3] < kClmin) {
1440 //printf("Exit zouts[1] = %d, zouts[3] = %d\n",zouts[1],zouts[3]);
1441 return;
1442 }
1443
1444 // Z distance bigger than pad - length
1445 if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) zouts[3] = 0;
1446
1447 Int_t breaktime = -1;
1448 Bool_t mbefore = kFALSE;
1449 Int_t cumul[kNtb][2];
1450 Int_t counts[2] = { 0, 0 };
1451
1452 if (zouts[3] >= 3) {
1453
1454 //
1455 // Find the break time allowing one chage on pad-rows
1456 // with maximal number of accepted clusters
1457 //
1458 fNChange = 1;
1459 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1460 cumul[i][0] = counts[0];
1461 cumul[i][1] = counts[1];
1462 if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++;
1463 if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++;
1464 }
1465 Int_t maxcount = 0;
1466 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1467 Int_t after = cumul[AliTRDtrackerV1::GetNTimeBins()][0] - cumul[i][0];
1468 Int_t before = cumul[i][1];
1469 if (after + before > maxcount) {
1470 maxcount = after + before;
1471 breaktime = i;
1472 mbefore = kFALSE;
1473 }
1474 after = cumul[AliTRDtrackerV1::GetNTimeBins()-1][1] - cumul[i][1];
1475 before = cumul[i][0];
1476 if (after + before > maxcount) {
1477 maxcount = after + before;
1478 breaktime = i;
1479 mbefore = kTRUE;
1480 }
1481 }
1482 breaktime -= 1;
1483 }
1484
1485 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1486 if (i > breaktime) allowedz[i] = mbefore ? zouts[2] : zouts[0];
1487 if (i <= breaktime) allowedz[i] = (!mbefore) ? zouts[2] : zouts[0];
1488 }
1489
1490 if (((allowedz[0] > allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] < 0)) ||
1491 ((allowedz[0] < allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] > 0))) {
1492 //
1493 // Tracklet z-direction not in correspondance with track z direction
1494 //
1495 fNChange = 0;
1496 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1497 allowedz[i] = zouts[0]; // Only longest taken
1498 }
1499 }
1500
1501 if (fNChange > 0) {
1502 //
1503 // Cross pad -row tracklet - take the step change into account
1504 //
1505 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1506 if (!fClusters[i]) continue;
1507 if(!fClusters[i]->IsInChamber()) continue;
1508 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
1509 // Residual y
dd8059a8 1510 //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] + GetTilt()*(fZ[i] - fZref[0]);
1511 yres[i] = fY[i] - GetTilt()*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
f29f13a6 1512// if (TMath::Abs(fZ[i] - fZProb) > 2) {
dd8059a8 1513// if (fZ[i] > fZProb) yres[i] += GetTilt() * GetPadLength();
1514// if (fZ[i] < fZProb) yres[i] -= GetTilt() * GetPadLength();
f29f13a6 1515 }
e3cf3d02 1516 }
1517 }
1518
1519 Double_t yres2[kNtb];
1520 Double_t mean;
1521 Double_t sigma;
1522 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1523 if (!fClusters[i]) continue;
1524 if(!fClusters[i]->IsInChamber()) continue;
1525 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
1526 yres2[fN2] = yres[i];
1527 fN2++;
1528 }
1529 if (fN2 < kClmin) {
1530 //printf("Exit fN2 < kClmin: fN2 = %d\n", fN2);
1531 fN2 = 0;
1532 return;
1533 }
1534 AliMathBase::EvaluateUni(fN2,yres2,mean,sigma, Int_t(fN2*kRatio-2.));
1535 if (sigma < sigmaexp * 0.8) {
1536 sigma = sigmaexp;
1537 }
1538 //Float_t fSigmaY = sigma;
1539
1540 // Reset sums
1541 sumw = 0;
1542 sumwx = 0;
1543 sumwx2 = 0;
1544 sumwy = 0;
1545 sumwxy = 0;
1546 sumwz = 0;
1547 sumwxz = 0;
1548
1549 fN2 = 0;
1550 Float_t fMeanz = 0;
1551 Float_t fMPads = 0;
1552 fUsable = 0;
1553 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1554 if (!fClusters[i]) continue;
1555 if (!fClusters[i]->IsInChamber()) continue;
1556 if (TMath::Abs(fZ[i] - allowedz[i]) > 2){fClusters[i] = 0x0; continue;}
1557 if (TMath::Abs(yres[i] - mean) > 4.0 * sigma){fClusters[i] = 0x0; continue;}
1558 SETBIT(fUsable,i);
1559 fN2++;
1560 fMPads += fClusters[i]->GetNPads();
1561 Float_t weight = 1.0;
1562 if (fClusters[i]->GetNPads() > 4) weight = 0.5;
1563 if (fClusters[i]->GetNPads() > 5) weight = 0.2;
1564
1565
1566 Double_t x = fX[i];
1567 //printf("x = %7.3f dy = %7.3f fit %7.3f\n", x, yres[i], fY[i]-yres[i]);
1568
1569 sumw += weight;
1570 sumwx += x * weight;
1571 sumwx2 += x*x * weight;
1572 sumwy += weight * yres[i];
1573 sumwxy += weight * (yres[i]) * x;
1574 sumwz += weight * fZ[i];
1575 sumwxz += weight * fZ[i] * x;
1576
1577 }
1578
1579 if (fN2 < kClmin){
1580 //printf("Exit fN2 < kClmin(2): fN2 = %d\n",fN2);
1581 fN2 = 0;
1582 return;
1583 }
1584 fMeanz = sumwz / sumw;
1585 Float_t correction = 0;
1586 if (fNChange > 0) {
1587 // Tracklet on boundary
1588 if (fMeanz < fZProb) correction = ycrosscor;
1589 if (fMeanz > fZProb) correction = -ycrosscor;
1590 }
1591
1592 Double_t det = sumw * sumwx2 - sumwx * sumwx;
1593 fYfit[0] = (sumwx2 * sumwy - sumwx * sumwxy) / det;
1594 fYfit[1] = (sumw * sumwxy - sumwx * sumwy) / det;
1595
1596 fS2Y = 0;
1597 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1598 if (!TESTBIT(fUsable,i)) continue;
1599 Float_t delta = yres[i] - fYfit[0] - fYfit[1] * fX[i];
1600 fS2Y += delta*delta;
1601 }
1602 fS2Y = TMath::Sqrt(fS2Y / Float_t(fN2-2));
1603 // TEMPORARY UNTIL covariance properly calculated
1604 fS2Y = TMath::Max(fS2Y, Float_t(.1));
1605
1606 fZfit[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det;
1607 fZfit[1] = (sumw * sumwxz - sumwx * sumwz) / det;
1608// fYfitR[0] += fYref[0] + correction;
1609// fYfitR[1] += fYref[1];
1610// fYfit[0] = fYfitR[0];
1611 fYfit[1] = -fYfit[1];
1612
1613 UpdateUsed();
f29f13a6 1614}*/
e3cf3d02 1615
e4f2f73d 1616//___________________________________________________________________
203967fc 1617void AliTRDseedV1::Print(Option_t *o) const
e4f2f73d 1618{
1619 //
1620 // Printing the seedstatus
1621 //
1622
b72f4eaf 1623 AliInfo(Form("Det[%3d] X0[%7.2f] Pad{L[%5.2f] W[%5.2f] Tilt[%+6.2f]}", fDet, fX0, GetPadLength(), GetPadWidth(), GetTilt()));
dd8059a8 1624 AliInfo(Form("N[%2d] Nused[%2d] Nshared[%2d] [%d]", GetN(), GetNUsed(), GetNShared(), fN));
b72f4eaf 1625 AliInfo(Form("FLAGS : RC[%c] Kink[%c] SA[%c]", IsRowCross()?'y':'n', IsKink()?'y':'n', IsStandAlone()?'y':'n'));
dd8059a8 1626
1627 Double_t cov[3], x=GetX();
1628 GetCovAt(x, cov);
1629 AliInfo(" | x[cm] | y[cm] | z[cm] | dydx | dzdx |");
1630 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 1631 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 1632
1633
1634 if(strcmp(o, "a")!=0) return;
1635
4dc4dc2e 1636 AliTRDcluster* const* jc = &fClusters[0];
8d2bec9e 1637 for(int ic=0; ic<kNclusters; ic++, jc++) {
4dc4dc2e 1638 if(!(*jc)) continue;
203967fc 1639 (*jc)->Print(o);
4dc4dc2e 1640 }
e4f2f73d 1641}
47d5d320 1642
203967fc 1643
1644//___________________________________________________________________
1645Bool_t AliTRDseedV1::IsEqual(const TObject *o) const
1646{
1647 // Checks if current instance of the class has the same essential members
1648 // as the given one
1649
1650 if(!o) return kFALSE;
1651 const AliTRDseedV1 *inTracklet = dynamic_cast<const AliTRDseedV1*>(o);
1652 if(!inTracklet) return kFALSE;
1653
1654 for (Int_t i = 0; i < 2; i++){
e3cf3d02 1655 if ( fYref[i] != inTracklet->fYref[i] ) return kFALSE;
1656 if ( fZref[i] != inTracklet->fZref[i] ) return kFALSE;
203967fc 1657 }
1658
e3cf3d02 1659 if ( fS2Y != inTracklet->fS2Y ) return kFALSE;
dd8059a8 1660 if ( GetTilt() != inTracklet->GetTilt() ) return kFALSE;
1661 if ( GetPadLength() != inTracklet->GetPadLength() ) return kFALSE;
203967fc 1662
8d2bec9e 1663 for (Int_t i = 0; i < kNclusters; i++){
e3cf3d02 1664// if ( fX[i] != inTracklet->GetX(i) ) return kFALSE;
1665// if ( fY[i] != inTracklet->GetY(i) ) return kFALSE;
1666// if ( fZ[i] != inTracklet->GetZ(i) ) return kFALSE;
1667 if ( fIndexes[i] != inTracklet->fIndexes[i] ) return kFALSE;
203967fc 1668 }
f29f13a6 1669// if ( fUsable != inTracklet->fUsable ) return kFALSE;
203967fc 1670
1671 for (Int_t i=0; i < 2; i++){
e3cf3d02 1672 if ( fYfit[i] != inTracklet->fYfit[i] ) return kFALSE;
1673 if ( fZfit[i] != inTracklet->fZfit[i] ) return kFALSE;
1674 if ( fLabels[i] != inTracklet->fLabels[i] ) return kFALSE;
203967fc 1675 }
1676
e3cf3d02 1677/* if ( fMeanz != inTracklet->GetMeanz() ) return kFALSE;
1678 if ( fZProb != inTracklet->GetZProb() ) return kFALSE;*/
3e778975 1679 if ( fN != inTracklet->fN ) return kFALSE;
1680 //if ( fNUsed != inTracklet->fNUsed ) return kFALSE;
e3cf3d02 1681 //if ( fFreq != inTracklet->GetFreq() ) return kFALSE;
1682 //if ( fNChange != inTracklet->GetNChange() ) return kFALSE;
203967fc 1683
e3cf3d02 1684 if ( fC != inTracklet->fC ) return kFALSE;
1685 //if ( fCC != inTracklet->GetCC() ) return kFALSE;
1686 if ( fChi2 != inTracklet->fChi2 ) return kFALSE;
203967fc 1687 // if ( fChi2Z != inTracklet->GetChi2Z() ) return kFALSE;
1688
e3cf3d02 1689 if ( fDet != inTracklet->fDet ) return kFALSE;
b25a5e9e 1690 if ( fPt != inTracklet->fPt ) return kFALSE;
e3cf3d02 1691 if ( fdX != inTracklet->fdX ) return kFALSE;
203967fc 1692
8d2bec9e 1693 for (Int_t iCluster = 0; iCluster < kNclusters; iCluster++){
203967fc 1694 AliTRDcluster *curCluster = fClusters[iCluster];
e3cf3d02 1695 AliTRDcluster *inCluster = inTracklet->fClusters[iCluster];
203967fc 1696 if (curCluster && inCluster){
1697 if (! curCluster->IsEqual(inCluster) ) {
1698 curCluster->Print();
1699 inCluster->Print();
1700 return kFALSE;
1701 }
1702 } else {
1703 // if one cluster exists, and corresponding
1704 // in other tracklet doesn't - return kFALSE
1705 if(curCluster || inCluster) return kFALSE;
1706 }
1707 }
1708 return kTRUE;
1709}
5d401b45 1710