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