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