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