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