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