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