]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliVertexerTracks.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Introducing a new AliVTrack class which will keep the features common for the ESD...
[u/mrichter/AliRoot.git] / STEER / AliVertexerTracks.cxx
CommitLineData
2d57349e 1/**************************************************************************
2 * Copyright(c) 1998-2003, 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 **************************************************************************/
15
16//-----------------------------------------------------------------
6d8df534 17// Implementation of the vertexer from tracks
2d57349e 18//
19// Origin: AliITSVertexerTracks
20// A.Dainese, Padova,
21// andrea.dainese@pd.infn.it
22// M.Masera, Torino,
23// massimo.masera@to.infn.it
24// Moved to STEER and adapted to ESD tracks:
6d8df534 25// F.Prino, Torino, prino@to.infn.it
2d57349e 26//-----------------------------------------------------------------
27
28//---- Root headers --------
a7ede274 29#include <TSystem.h>
30#include <TDirectory.h>
31#include <TFile.h>
2d57349e 32#include <TTree.h>
ec1be5d5 33#include <TMatrixD.h>
2d57349e 34//---- AliRoot headers -----
919e537f 35#include "AliLog.h"
2d57349e 36#include "AliStrLine.h"
6d8df534 37#include "AliExternalTrackParam.h"
af885e0f 38#include "AliESDEvent.h"
2d57349e 39#include "AliESDtrack.h"
6d8df534 40#include "AliVertexerTracks.h"
2d57349e 41
42ClassImp(AliVertexerTracks)
43
44
45//----------------------------------------------------------------------------
46AliVertexerTracks::AliVertexerTracks():
07680cae 47TObject(),
48fVert(),
49fCurrentVertex(0),
f09c879d 50fMode(0),
dc3719f3 51fFieldkG(-999.),
6d8df534 52fTrkArraySel(),
53fIdSel(0),
07680cae 54fTrksToSkip(0),
55fNTrksToSkip(0),
f09c879d 56fConstraint(kFALSE),
57fOnlyFitter(kFALSE),
58fMinTracks(1),
6b29d399 59fMinClusters(5),
f09c879d 60fDCAcut(0.1),
61fDCAcutIter0(0.1),
62fNSigma(3.),
fee9ef0d 63fMaxd0z0(0.5),
f09c879d 64fMinDetFitter(100.),
65fMaxTgl(1000.),
50ff0bcd 66fITSrefit(kTRUE),
6b29d399 67fFiducialR(3.),
68fFiducialZ(30.),
3b768a1d 69fnSigmaForUi00(1.5),
f09c879d 70fAlgo(1)
2d57349e 71{
72//
73// Default constructor
74//
3b768a1d 75 SetVtxStart();
07680cae 76 SetVtxStartSigma();
2d57349e 77}
dc3719f3 78//----------------------------------------------------------------------------
79AliVertexerTracks::AliVertexerTracks(Double_t fieldkG):
07680cae 80TObject(),
81fVert(),
82fCurrentVertex(0),
f09c879d 83fMode(0),
84fFieldkG(fieldkG),
6d8df534 85fTrkArraySel(),
86fIdSel(0),
07680cae 87fTrksToSkip(0),
88fNTrksToSkip(0),
f09c879d 89fConstraint(kFALSE),
90fOnlyFitter(kFALSE),
91fMinTracks(1),
6b29d399 92fMinClusters(5),
f09c879d 93fDCAcut(0.1),
94fDCAcutIter0(0.1),
95fNSigma(3.),
fee9ef0d 96fMaxd0z0(0.5),
f09c879d 97fMinDetFitter(100.),
98fMaxTgl(1000.),
50ff0bcd 99fITSrefit(kTRUE),
6b29d399 100fFiducialR(3.),
101fFiducialZ(30.),
3b768a1d 102fnSigmaForUi00(1.5),
f09c879d 103fAlgo(1)
2d57349e 104{
105//
106// Standard constructor
107//
3b768a1d 108 SetVtxStart();
07680cae 109 SetVtxStartSigma();
6b29d399 110 SetTPCMode();
2d57349e 111}
112//-----------------------------------------------------------------------------
fee9ef0d 113AliVertexerTracks::~AliVertexerTracks()
114{
2d57349e 115 // Default Destructor
eae1677e 116 // The objects pointed by the following pointer are not owned
2d57349e 117 // by this class and are not deleted
ec1be5d5 118 fCurrentVertex = 0;
146c29df 119 if(fTrksToSkip) { delete [] fTrksToSkip; fTrksToSkip=NULL; }
6d8df534 120 if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
ec1be5d5 121}
c5e3e5d1 122//----------------------------------------------------------------------------
af885e0f 123AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESDEvent *esdEvent)
ec1be5d5 124{
c5e3e5d1 125//
ec1be5d5 126// Primary vertex for current ESD event
eae1677e 127// (Two iterations:
fee9ef0d 128// 1st with 5*fNSigma*sigma cut w.r.t. to initial vertex
129// + cut on sqrt(d0d0+z0z0) if fConstraint=kFALSE
130// 2nd with fNSigma*sigma cut w.r.t. to vertex found in 1st iteration)
07680cae 131//
132 fCurrentVertex = 0;
133
dc3719f3 134 // accept 1-track case only if constraint is available
135 if(!fConstraint && fMinTracks==1) fMinTracks=2;
136
fee9ef0d 137 // read tracks from ESD
6d8df534 138 Int_t nTrks = (Int_t)esdEvent->GetNumberOfTracks();
f09c879d 139 if(nTrks<fMinTracks) {
6d8df534 140 TooFewTracks();
80d36431 141 return fCurrentVertex;
142 }
6d8df534 143 //
713e5fbf 144
a7ede274 145 TDirectory * olddir = gDirectory;
f09c879d 146 TFile *f = 0;
147 if(nTrks>500) f = new TFile("VertexerTracks.root","recreate");
6d8df534 148 TObjArray trkArrayOrig(nTrks);
149 UShort_t *idOrig = new UShort_t[nTrks];
07680cae 150
6d8df534 151 Int_t nTrksOrig=0;
f09c879d 152 AliExternalTrackParam *t=0;
153 // loop on ESD tracks
6d8df534 154 for(Int_t i=0; i<nTrks; i++) {
155 AliESDtrack *esdt = esdEvent->GetTrack(i);
eae1677e 156 // check tracks to skip
6d8df534 157 Bool_t skipThis = kFALSE;
eae1677e 158 for(Int_t j=0; j<fNTrksToSkip; j++) {
6d8df534 159 if(esdt->GetID()==fTrksToSkip[j]) {
919e537f 160 AliDebug(1,Form("skipping track: %d",i));
eae1677e 161 skipThis = kTRUE;
162 }
163 }
a7ede274 164 if(skipThis) continue;
f09c879d 165
6b29d399 166 // check number of clusters in ITS or TPC
167 if(esdt->GetNcls(fMode) < fMinClusters) continue;
168
f09c879d 169 if(fMode==0) { // ITS mode
170 if(fITSrefit && !(esdt->GetStatus()&AliESDtrack::kITSrefit)) continue;
f09c879d 171 Double_t x,p[5],cov[15];
172 esdt->GetExternalParameters(x,p);
173 esdt->GetExternalCovariance(cov);
174 t = new AliExternalTrackParam(x,esdt->GetAlpha(),p,cov);
175 } else if(fMode==1) { // TPC mode
176 t = (AliExternalTrackParam*)esdt->GetTPCInnerParam();
177 if(!t) continue;
178 Double_t radius = 2.8; //something less than the beam pipe radius
179 if(!PropagateTrackTo(t,radius)) continue;
50ff0bcd 180 }
6d8df534 181 trkArrayOrig.AddLast(t);
182 idOrig[nTrksOrig]=(UShort_t)esdt->GetID();
183 nTrksOrig++;
f09c879d 184 } // end loop on ESD tracks
fee9ef0d 185
f09c879d 186 // call method that will reconstruct the vertex
6d8df534 187 FindPrimaryVertex(&trkArrayOrig,idOrig);
188
f09c879d 189 if(fMode==0) trkArrayOrig.Delete();
6d8df534 190 delete [] idOrig; idOrig=NULL;
f09c879d 191
192 if(f) {
193 f->Close(); delete f; f = NULL;
194 gSystem->Unlink("VertexerTracks.root");
195 olddir->cd();
196 }
6d8df534 197
198 return fCurrentVertex;
199}
200//----------------------------------------------------------------------------
201AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
202 UShort_t *idOrig)
203{
204//
205// Primary vertex using the AliExternalTrackParam's in the TObjArray.
206// idOrig must contain the track IDs from AliESDtrack::GetID()
207// (Two iterations:
208// 1st with 5*fNSigma*sigma cut w.r.t. to initial vertex
209// + cut on sqrt(d0d0+z0z0) if fConstraint=kFALSE
210// 2nd with fNSigma*sigma cut w.r.t. to vertex found in 1st iteration)
211//
212 fCurrentVertex = 0;
213
214 // accept 1-track case only if constraint is available
215 if(!fConstraint && fMinTracks==1) fMinTracks=2;
216
f09c879d 217 // read tracks from array
6d8df534 218 Int_t nTrksOrig = (Int_t)trkArrayOrig->GetEntriesFast();
919e537f 219 AliDebug(1,Form("Initial number of tracks: %d",nTrksOrig));
f09c879d 220 if(nTrksOrig<fMinTracks) {
919e537f 221 AliDebug(1,"TooFewTracks");
6d8df534 222 TooFewTracks();
223 return fCurrentVertex;
224 }
225
fee9ef0d 226 // If fConstraint=kFALSE
227 // run VertexFinder(1) to get rough estimate of initVertex (x,y)
228 if(!fConstraint) {
6d8df534 229 // fill fTrkArraySel, for VertexFinder()
230 fIdSel = new UShort_t[nTrksOrig];
231 PrepareTracks(*trkArrayOrig,idOrig,0);
232 if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
233 Double_t cutsave = fDCAcut;
f09c879d 234 fDCAcut = fDCAcutIter0;
235 VertexFinder(1); // using weights, cutting dca < fDCAcutIter0
fee9ef0d 236 fDCAcut = cutsave;
fee9ef0d 237 if(fVert.GetNContributors()>0) {
238 fVert.GetXYZ(fNominalPos);
239 fNominalPos[0] = fVert.GetXv();
240 fNominalPos[1] = fVert.GetYv();
241 fNominalPos[2] = fVert.GetZv();
919e537f 242 AliDebug(1,Form("No mean vertex: VertexFinder gives (%f, %f, %f)",fNominalPos[0],fNominalPos[1],fNominalPos[2]));
fee9ef0d 243 } else {
244 fNominalPos[0] = 0.;
245 fNominalPos[1] = 0.;
246 fNominalPos[2] = 0.;
919e537f 247 AliDebug(1,"No mean vertex and VertexFinder failed");
0c69be49 248 }
07680cae 249 }
fee9ef0d 250
251 // TWO ITERATIONS:
252 //
253 // ITERATION 1
254 // propagate tracks to fNominalPos vertex
255 // preselect them:
256 // if(constraint) reject for |d0|>5*fNSigma*sigma w.r.t. fNominal... vertex
257 // else reject for |d0|\oplus|z0| > 5 mm w.r.t. fNominal... vertex
07680cae 258 // ITERATION 2
259 // propagate tracks to best between initVertex and fCurrentVertex
260 // preselect tracks (reject for |d0|>fNSigma*sigma w.r.t. best
261 // between initVertex and fCurrentVertex)
6d8df534 262 for(Int_t iter=1; iter<=2; iter++) {
263 if(fOnlyFitter && iter==1) continue;
264 if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
265 fIdSel = new UShort_t[nTrksOrig];
266 Int_t nTrksSel = PrepareTracks(*trkArrayOrig,idOrig,iter);
919e537f 267 AliDebug(1,Form("N tracks selected in iteration %d: %d",iter,nTrksSel));
6d8df534 268 if(nTrksSel < fMinTracks) {
269 TooFewTracks();
fee9ef0d 270 return fCurrentVertex;
271 }
07680cae 272
fee9ef0d 273 // vertex finder
274 if(!fOnlyFitter) {
6d8df534 275 if(nTrksSel==1) {
919e537f 276 AliDebug(1,"Just one track");
fee9ef0d 277 OneTrackVertFinder();
6d8df534 278 } else {
fee9ef0d 279 switch (fAlgo) {
280 case 1: StrLinVertexFinderMinDist(1); break;
281 case 2: StrLinVertexFinderMinDist(0); break;
282 case 3: HelixVertexFinder(); break;
283 case 4: VertexFinder(1); break;
284 case 5: VertexFinder(0); break;
919e537f 285 default: printf("Wrong algorithm"); break;
fee9ef0d 286 }
287 }
919e537f 288 AliDebug(1," Vertex finding completed");
0c69be49 289 }
07680cae 290
fee9ef0d 291 // vertex fitter
6d8df534 292 VertexFitter();
fee9ef0d 293 } // end loop on the two iterations
07680cae 294
fee9ef0d 295
50ff0bcd 296 // set indices of used tracks
297 UShort_t *indices = 0;
50ff0bcd 298 if(fCurrentVertex->GetNContributors()>0) {
4d023a8f 299 Int_t nIndices = (Int_t)fTrkArraySel.GetEntriesFast();
300 indices = new UShort_t[nIndices];
301 for(Int_t jj=0; jj<nIndices; jj++)
6d8df534 302 indices[jj] = fIdSel[jj];
4d023a8f 303 fCurrentVertex->SetIndices(nIndices,indices);
50ff0bcd 304 }
6d8df534 305 delete [] indices; indices=NULL;
306 //
fee9ef0d 307
6d8df534 308 // set vertex title
309 TString title="VertexerTracksNoConstraint";
310 if(fConstraint) {
311 title="VertexerTracksWithConstraint";
312 if(fOnlyFitter) title.Append("OnlyFitter");
313 }
314 fCurrentVertex->SetTitle(title.Data());
315 //
a7ede274 316
919e537f 317 AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
50ff0bcd 318
6d8df534 319 // clean up
320 delete [] fIdSel; fIdSel=NULL;
321 fTrkArraySel.Delete();
146c29df 322 if(fTrksToSkip) { delete [] fTrksToSkip; fTrksToSkip=NULL; }
6d8df534 323 //
fee9ef0d 324
07680cae 325 return fCurrentVertex;
326}
50ff0bcd 327//------------------------------------------------------------------------
fee9ef0d 328Double_t AliVertexerTracks::GetDeterminant3X3(Double_t matr[][3])
329{
50ff0bcd 330 //
331 Double_t det=matr[0][0]*matr[1][1]*matr[2][2]-matr[0][0]*matr[1][2]*matr[2][1]-matr[0][1]*matr[1][0]*matr[2][2]+matr[0][1]*matr[1][2]*matr[2][0]+matr[0][2]*matr[1][0]*matr[2][1]-matr[0][2]*matr[1][1]*matr[2][0];
332 return det;
333}
334//-------------------------------------------------------------------------
fee9ef0d 335void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t (*m)[3],Double_t *d)
336{
50ff0bcd 337 //
338 Double_t x12=p0[0]-p1[0];
339 Double_t y12=p0[1]-p1[1];
340 Double_t z12=p0[2]-p1[2];
341 Double_t kk=x12*x12+y12*y12+z12*z12;
342 m[0][0]=2-2/kk*x12*x12;
343 m[0][1]=-2/kk*x12*y12;
344 m[0][2]=-2/kk*x12*z12;
345 m[1][0]=-2/kk*x12*y12;
346 m[1][1]=2-2/kk*y12*y12;
347 m[1][2]=-2/kk*y12*z12;
348 m[2][0]=-2/kk*x12*z12;
349 m[2][1]=-2*y12*z12;
350 m[2][2]=2-2/kk*z12*z12;
351 d[0]=2*p0[0]-2/kk*p0[0]*x12*x12-2/kk*p0[2]*x12*z12-2/kk*p0[1]*x12*y12;
352 d[1]=2*p0[1]-2/kk*p0[1]*y12*y12-2/kk*p0[0]*x12*y12-2/kk*p0[2]*z12*y12;
353 d[2]=2*p0[2]-2/kk*p0[2]*z12*z12-2/kk*p0[0]*x12*z12-2/kk*p0[1]*z12*y12;
c5e3e5d1 354
50ff0bcd 355}
356//--------------------------------------------------------------------------
fee9ef0d 357void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t *sigmasq,Double_t (*m)[3],Double_t *d)
358{
50ff0bcd 359 //
360 Double_t x12=p1[0]-p0[0];
361 Double_t y12=p1[1]-p0[1];
362 Double_t z12=p1[2]-p0[2];
c5e3e5d1 363
50ff0bcd 364 Double_t den= x12*x12*sigmasq[1]*sigmasq[2]+y12*y12*sigmasq[0]*sigmasq[2]+z12*z12*sigmasq[0]*sigmasq[1];
c5e3e5d1 365
50ff0bcd 366 Double_t kk= 2*(x12*x12/sigmasq[0]+y12*y12/sigmasq[1]+z12*z12/sigmasq[2]);
c5e3e5d1 367
50ff0bcd 368 Double_t cc[3];
369 cc[0]=-x12/sigmasq[0];
370 cc[1]=-y12/sigmasq[1];
371 cc[2]=-z12/sigmasq[2];
ec1be5d5 372
50ff0bcd 373 Double_t ww=(-p0[0]*x12*sigmasq[1]*sigmasq[2]-p0[1]*y12*sigmasq[0]*sigmasq[2]-p0[2]*z12*sigmasq[0]*sigmasq[1])/den;
ec1be5d5 374
50ff0bcd 375 Double_t ss= -p0[0]*cc[0]-p0[1]*cc[1]-p0[2]*cc[2];
2d57349e 376
50ff0bcd 377 Double_t aa[3];
378 aa[0]=x12*sigmasq[1]*sigmasq[2]/den;
379 aa[1]=y12*sigmasq[0]*sigmasq[2]/den;
380 aa[2]=z12*sigmasq[0]*sigmasq[1]/den;
2d57349e 381
50ff0bcd 382 m[0][0]=aa[0]*(aa[0]*kk+2*cc[0])+2*cc[0]*aa[0]+2/sigmasq[0];
383 m[0][1]=aa[1]*(aa[0]*kk+2*cc[0])+2*cc[1]*aa[0];
384 m[0][2]=aa[2]*(aa[0]*kk+2*cc[0])+2*cc[2]*aa[0];
2d57349e 385
50ff0bcd 386 m[1][0]=aa[0]*(aa[1]*kk+2*cc[1])+2*cc[0]*aa[1];
387 m[1][1]=aa[1]*(aa[1]*kk+2*cc[1])+2*cc[1]*aa[1]+2/sigmasq[1];
388 m[1][2]=aa[2]*(aa[1]*kk+2*cc[1])+2*cc[2]*aa[1];
2d57349e 389
50ff0bcd 390 m[2][0]=aa[0]*(aa[2]*kk+2*cc[2])+2*cc[0]*aa[2];
391 m[2][1]=aa[1]*(aa[2]*kk+2*cc[2])+2*cc[1]*aa[2];
392 m[2][2]=aa[2]*(aa[2]*kk+2*cc[2])+2*cc[2]*aa[2]+2/sigmasq[2];
08df6187 393
50ff0bcd 394 d[0]=-ww*(aa[0]*kk+2*cc[0])-2*ss*aa[0]+2*p0[0]/sigmasq[0];
395 d[1]=-ww*(aa[1]*kk+2*cc[1])-2*ss*aa[1]+2*p0[1]/sigmasq[1];
396 d[2]=-ww*(aa[2]*kk+2*cc[2])-2*ss*aa[2]+2*p0[2]/sigmasq[2];
2d57349e 397
2d57349e 398 }
50ff0bcd 399//--------------------------------------------------------------------------
fee9ef0d 400Double_t AliVertexerTracks::GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t *x0)
401{
50ff0bcd 402 //
403 Double_t x12=p0[0]-p1[0];
404 Double_t y12=p0[1]-p1[1];
405 Double_t z12=p0[2]-p1[2];
406 Double_t x10=p0[0]-x0[0];
407 Double_t y10=p0[1]-x0[1];
408 Double_t z10=p0[2]-x0[2];
3446476f 409 // return ((x10*x10+y10*y10+z10*z10)*(x12*x12+y12*y12+z12*z12)-(x10*x12+y10*y12+z10*z12)*(x10*x12+y10*y12+z10*z12))/(x12*x12+y12*y12+z12*z12);
3446476f 410 return ((y10*z12-z10*y12)*(y10*z12-z10*y12)+
411 (z10*x12-x10*z12)*(z10*x12-x10*z12)+
412 (x10*y12-y10*x12)*(x10*y12-y10*x12))
413 /(x12*x12+y12*y12+z12*z12);
2d57349e 414}
415//---------------------------------------------------------------------------
fee9ef0d 416void AliVertexerTracks::OneTrackVertFinder()
417{
0c69be49 418 // find vertex for events with 1 track, using DCA to nominal beam axis
919e537f 419 AliDebug(1,Form("Number of prepared tracks =%d - Call OneTrackVertFinder",fTrkArraySel.GetEntries()));
6d8df534 420 AliExternalTrackParam *track1;
421 track1 = (AliExternalTrackParam*)fTrkArraySel.At(0);
0c69be49 422 Double_t alpha=track1->GetAlpha();
423 Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
424 Double_t pos[3],dir[3];
f09c879d 425 track1->GetXYZAt(mindist,GetFieldkG(),pos);
426 track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
0c69be49 427 AliStrLine *line1 = new AliStrLine(pos,dir);
428 Double_t p1[3]={fNominalPos[0],fNominalPos[1],0.};
429 Double_t p2[3]={fNominalPos[0],fNominalPos[1],10.};
430 AliStrLine *zeta=new AliStrLine(p1,p2,kTRUE);
431 Double_t crosspoint[3]={0.,0.,0.};
432 Double_t sigma=999.;
433 Int_t nContrib=-1;
434 Int_t retcode = zeta->Cross(line1,crosspoint);
435 if(retcode>=0){
436 sigma=line1->GetDistFromPoint(crosspoint);
437 nContrib=1;
438 }
439 delete zeta;
440 delete line1;
441 fVert.SetXYZ(crosspoint);
442 fVert.SetDispersion(sigma);
443 fVert.SetNContributors(nContrib);
444}
445//---------------------------------------------------------------------------
fee9ef0d 446void AliVertexerTracks::HelixVertexFinder()
447{
2d57349e 448 // Get estimate of vertex position in (x,y) from tracks DCA
449
450
451 Double_t initPos[3];
452 initPos[2] = 0.;
453 for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
2d57349e 454
6d8df534 455 Int_t nacc = (Int_t)fTrkArraySel.GetEntriesFast();
2d57349e 456
457 Double_t aver[3]={0.,0.,0.};
458 Double_t averquad[3]={0.,0.,0.};
459 Double_t sigmaquad[3]={0.,0.,0.};
460 Double_t sigma=0;
461 Int_t ncombi = 0;
6d8df534 462 AliExternalTrackParam *track1;
463 AliExternalTrackParam *track2;
2d57349e 464 Double_t distCA;
6d8df534 465 Double_t x;
2d57349e 466 Double_t alpha, cs, sn;
467 Double_t crosspoint[3];
468 for(Int_t i=0; i<nacc; i++){
6d8df534 469 track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
2d57349e 470
471
472 for(Int_t j=i+1; j<nacc; j++){
6d8df534 473 track2 = (AliExternalTrackParam*)fTrkArraySel.At(j);
2d57349e 474
f09c879d 475 distCA=track2->PropagateToDCA(track1,GetFieldkG());
2d57349e 476 if(fDCAcut<=0 ||(fDCAcut>0&&distCA<fDCAcut)){
6d8df534 477 x=track1->GetX();
2d57349e 478 alpha=track1->GetAlpha();
479 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
6d8df534 480 Double_t x1=x*cs - track1->GetY()*sn;
481 Double_t y1=x*sn + track1->GetY()*cs;
482 Double_t z1=track1->GetZ();
483
2d57349e 484 Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
6d8df534 485 x=track2->GetX();
2d57349e 486 alpha=track2->GetAlpha();
487 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
6d8df534 488 Double_t x2=x*cs - track2->GetY()*sn;
489 Double_t y2=x*sn + track2->GetY()*cs;
490 Double_t z2=track2->GetZ();
2d57349e 491 Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
492 Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
493 Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
494 Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
495 Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
496 crosspoint[0]=wx1*x1 + wx2*x2;
497 crosspoint[1]=wy1*y1 + wy2*y2;
498 crosspoint[2]=wz1*z1 + wz2*z2;
499
500 ncombi++;
501 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
502 for(Int_t jj=0;jj<3;jj++)averquad[jj]+=(crosspoint[jj]*crosspoint[jj]);
503 }
504 }
505
506 }
507 if(ncombi>0){
508 for(Int_t jj=0;jj<3;jj++){
509 initPos[jj] = aver[jj]/ncombi;
510 averquad[jj]/=ncombi;
511 sigmaquad[jj]=averquad[jj]-initPos[jj]*initPos[jj];
512 sigma+=sigmaquad[jj];
513 }
514 sigma=TMath::Sqrt(TMath::Abs(sigma));
515 }
516 else {
517 Warning("HelixVertexFinder","Finder did not succed");
518 sigma=999;
519 }
520 fVert.SetXYZ(initPos);
521 fVert.SetDispersion(sigma);
522 fVert.SetNContributors(ncombi);
523}
50ff0bcd 524//----------------------------------------------------------------------------
6d8df534 525Int_t AliVertexerTracks::PrepareTracks(TObjArray &trkArrayOrig,
526 UShort_t *idOrig,
527 Int_t optImpParCut)
fee9ef0d 528{
50ff0bcd 529//
530// Propagate tracks to initial vertex position and store them in a TObjArray
531//
919e537f 532 AliDebug(1," PrepareTracks()");
6d8df534 533
534 Int_t nTrksOrig = (Int_t)trkArrayOrig.GetEntriesFast();
535 Int_t nTrksSel = 0;
fee9ef0d 536 Double_t maxd0rphi;
50ff0bcd 537 Double_t sigmaCurr[3];
fee9ef0d 538 Double_t normdistx,normdisty;
6d8df534 539 Double_t d0z0[2],covd0z0[3];
43c9dae1 540 Double_t sigmad0;
50ff0bcd 541
542 AliESDVertex *initVertex = new AliESDVertex(fNominalPos,fNominalCov,1,1);
543
6d8df534 544 if(!fTrkArraySel.IsEmpty()) fTrkArraySel.Delete();
50ff0bcd 545
6d8df534 546 AliExternalTrackParam *track=0;
50ff0bcd 547
f09c879d 548 // loop on tracks
6d8df534 549 for(Int_t i=0; i<nTrksOrig; i++) {
550 track = new AliExternalTrackParam(*(AliExternalTrackParam*)trkArrayOrig.At(i));
f09c879d 551 // tgl cut
552 if(TMath::Abs(track->GetTgl())>fMaxTgl) {
919e537f 553 AliDebug(1,Form(" rejecting track with tgl = %f",track->GetTgl()));
f09c879d 554 delete track; continue;
6d8df534 555 }
50ff0bcd 556
18eba197 557 Bool_t propagateOK = kFALSE;
50ff0bcd 558 // propagate track to vertex
43c9dae1 559 if(optImpParCut<2 || fOnlyFitter) { // optImpParCut==1 or 0
18eba197 560 propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
fee9ef0d 561 } else { // optImpParCut==2
50ff0bcd 562 fCurrentVertex->GetSigmaXYZ(sigmaCurr);
fee9ef0d 563 normdistx = TMath::Abs(fCurrentVertex->GetXv()-fNominalPos[0])/TMath::Sqrt(sigmaCurr[0]*sigmaCurr[0]+fNominalCov[0]);
80d36431 564 normdisty = TMath::Abs(fCurrentVertex->GetYv()-fNominalPos[1])/TMath::Sqrt(sigmaCurr[1]*sigmaCurr[1]+fNominalCov[2]);
fee9ef0d 565 if(normdistx < 3. && normdisty < 3. &&
80d36431 566 (sigmaCurr[0]+sigmaCurr[1])<(TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2]))) {
18eba197 567 propagateOK = track->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
50ff0bcd 568 } else {
18eba197 569 propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
50ff0bcd 570 }
571 }
572
18eba197 573 if(!propagateOK) {
919e537f 574 AliDebug(1," rejected");
18eba197 575 delete track; continue;
576 }
577
43c9dae1 578 sigmad0 = TMath::Sqrt(covd0z0[0]);
579 maxd0rphi = fNSigma*sigmad0;
fee9ef0d 580 if(optImpParCut==1) maxd0rphi *= 5.;
6b29d399 581 maxd0rphi = TMath::Min(maxd0rphi,fFiducialR);
43c9dae1 582 //sigmad0z0 = TMath::Sqrt(covd0z0[0]+covd0z0[2]); // for future improvement
583 //maxd0z0 = 10.*fNSigma*sigmad0z0;
fee9ef0d 584
919e537f 585 AliDebug(1,Form("trk %d; id %d; |d0| = %f; d0 cut = %f; |z0| = %f; |d0|oplus|z0| = %f; d0z0 cut = %f",i,(Int_t)idOrig[i],TMath::Abs(d0z0[0]),maxd0rphi,TMath::Abs(d0z0[1]),TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1]),fMaxd0z0));
f09c879d 586
587
588 //---- track selection based on impact parameters ----//
589
590 // always reject tracks outside fiducial volume
6b29d399 591 if(TMath::Abs(d0z0[0])>fFiducialR || TMath::Abs(d0z0[1])>fFiducialZ) {
919e537f 592 AliDebug(1," rejected");
f09c879d 593 delete track; continue;
594 }
595
596 // during iterations 1 and 2 , reject tracks with d0rphi > maxd0rphi
597 if(optImpParCut>0 && TMath::Abs(d0z0[0])>maxd0rphi) {
919e537f 598 AliDebug(1," rejected");
f09c879d 599 delete track; continue;
600 }
fee9ef0d 601
43c9dae1 602 // if fConstraint=kFALSE, during iterations 1 and 2 ||
603 // if fConstraint=kTRUE, during iteration 2,
f09c879d 604 // select tracks with d0oplusz0 < fMaxd0z0
43c9dae1 605 if((!fConstraint && optImpParCut>0 && fVert.GetNContributors()>0) ||
606 ( fConstraint && optImpParCut==2)) {
607 if(nTrksOrig>=3 &&
f09c879d 608 TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1])>fMaxd0z0) {
919e537f 609 AliDebug(1," rejected");
fee9ef0d 610 delete track; continue;
611 }
612 }
43c9dae1 613
f09c879d 614 // track passed all selections
6d8df534 615 fTrkArraySel.AddLast(track);
616 fIdSel[nTrksSel] = idOrig[i];
617 nTrksSel++;
f09c879d 618 } // end loop on tracks
50ff0bcd 619
620 delete initVertex;
621
6d8df534 622 return nTrksSel;
50ff0bcd 623}
f09c879d 624//----------------------------------------------------------------------------
625Bool_t AliVertexerTracks::PropagateTrackTo(AliExternalTrackParam *track,
626 Double_t xToGo) {
627 //----------------------------------------------------------------
628 // COPIED from AliTracker
629 //
630 // Propagates the track to the plane X=xk (cm).
631 //
632 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
633 //----------------------------------------------------------------
634
635 const Double_t kEpsilon = 0.00001;
636 Double_t xpos = track->GetX();
637 Double_t dir = (xpos<xToGo) ? 1. : -1.;
638 Double_t maxStep = 5;
639 Double_t maxSnp = 0.8;
640 //
641 while ( (xToGo-xpos)*dir > kEpsilon){
642 Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
643 Double_t x = xpos+step;
644 Double_t xyz0[3],xyz1[3];
645 track->GetXYZ(xyz0); //starting global position
646
647 if(!track->GetXYZAt(x,GetFieldkG(),xyz1)) return kFALSE; // no prolongation
648 xyz1[2]+=kEpsilon; // waiting for bug correction in geo
649
650 if(TMath::Abs(track->GetSnpAt(x,GetFieldkG())) >= maxSnp) return kFALSE;
651 if(!track->PropagateTo(x,GetFieldkG())) return kFALSE;
652
653 if(TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
654 track->GetXYZ(xyz0); // global position
655 Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]);
656 //
657 Double_t ca=TMath::Cos(alphan-track->GetAlpha()),
658 sa=TMath::Sin(alphan-track->GetAlpha());
659 Double_t sf=track->GetSnp(), cf=TMath::Sqrt(1.- sf*sf);
660 Double_t sinNew = sf*ca - cf*sa;
661 if(TMath::Abs(sinNew) >= maxSnp) return kFALSE;
662 if(!track->Rotate(alphan)) return kFALSE;
663
664 xpos = track->GetX();
665 }
666 return kTRUE;
667}
50ff0bcd 668//---------------------------------------------------------------------------
fee9ef0d 669AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
6d8df534 670 TObjArray *trkArray,
671 UShort_t *id,
fee9ef0d 672 Float_t *diamondxy)
673{
50ff0bcd 674//
fee9ef0d 675// Removes tracks in trksTree from fit of inVtx
50ff0bcd 676//
fee9ef0d 677
146c29df 678 if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraint")) {
919e537f 679 printf("ERROR: primary vertex has no constraint: cannot remove tracks");
146c29df 680 return 0x0;
681 }
fee9ef0d 682 if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
919e537f 683 printf("WARNING: result of tracks' removal will be only approximately correct");
fee9ef0d 684
685 TMatrixD rv(3,1);
686 rv(0,0) = inVtx->GetXv();
687 rv(1,0) = inVtx->GetYv();
688 rv(2,0) = inVtx->GetZv();
689 TMatrixD vV(3,3);
690 Double_t cov[6];
691 inVtx->GetCovMatrix(cov);
692 vV(0,0) = cov[0];
693 vV(0,1) = cov[1]; vV(1,0) = cov[1];
694 vV(1,1) = cov[2];
695 vV(0,2) = cov[3]; vV(2,0) = cov[3];
696 vV(1,2) = cov[4]; vV(2,1) = cov[4];
697 vV(2,2) = cov[5];
698
699 TMatrixD sumWi(TMatrixD::kInverted,vV);
700 TMatrixD sumWiri(sumWi,TMatrixD::kMult,rv);
701
dad616ce 702 Int_t nUsedTrks = inVtx->GetNIndices();
fee9ef0d 703 Double_t chi2 = inVtx->GetChi2();
704
6d8df534 705 AliExternalTrackParam *track = 0;
706 Int_t ntrks = (Int_t)trkArray->GetEntriesFast();
fee9ef0d 707 for(Int_t i=0;i<ntrks;i++) {
6d8df534 708 track = (AliExternalTrackParam*)trkArray->At(i);
709 if(!inVtx->UsesTrack(id[i])) {
919e537f 710 printf("track %d was not used in vertex fit",id[i]);
fee9ef0d 711 continue;
712 }
713 Double_t alpha = track->GetAlpha();
714 Double_t xl = diamondxy[0]*TMath::Cos(alpha)+diamondxy[1]*TMath::Sin(alpha);
6d8df534 715 track->PropagateTo(xl,GetFieldkG());
fee9ef0d 716 // vector of track global coordinates
717 TMatrixD ri(3,1);
718 // covariance matrix of ri
719 TMatrixD wWi(3,3);
720
721 // get space point from track
722 if(!TrackToPoint(track,ri,wWi)) continue;
723
724 TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
725
726 sumWi -= wWi;
727 sumWiri -= wWiri;
728
729 // track chi2
730 TMatrixD deltar = rv; deltar -= ri;
731 TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
732 Double_t chi2i = deltar(0,0)*wWideltar(0,0)+
733 deltar(1,0)*wWideltar(1,0)+
734 deltar(2,0)*wWideltar(2,0);
735 // remove from total chi2
736 chi2 -= chi2i;
737
738 nUsedTrks--;
739 if(nUsedTrks<2) {
919e537f 740 printf("Trying to remove too many tracks!");
fee9ef0d 741 return 0x0;
742 }
743 }
744
745 TMatrixD rvnew(3,1);
746 TMatrixD vVnew(3,3);
747
748 // new inverted of weights matrix
749 TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
750 vVnew = invsumWi;
751 // new position of primary vertex
752 rvnew.Mult(vVnew,sumWiri);
753
754 Double_t position[3];
755 position[0] = rvnew(0,0);
756 position[1] = rvnew(1,0);
757 position[2] = rvnew(2,0);
758 cov[0] = vVnew(0,0);
759 cov[1] = vVnew(0,1);
760 cov[2] = vVnew(1,1);
761 cov[3] = vVnew(0,2);
762 cov[4] = vVnew(1,2);
763 cov[5] = vVnew(2,2);
764
765 // store data in the vertex object
766 AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,nUsedTrks);
767 outVtx->SetTitle(inVtx->GetTitle());
fee9ef0d 768 UShort_t *inindices = inVtx->GetIndices();
dad616ce 769 Int_t nIndices = nUsedTrks;
4d023a8f 770 UShort_t *outindices = new UShort_t[nIndices];
fee9ef0d 771 Int_t j=0;
fee9ef0d 772 for(Int_t k=0; k<inVtx->GetNIndices(); k++) {
6d8df534 773 Bool_t copyindex=kTRUE;
fee9ef0d 774 for(Int_t l=0; l<ntrks; l++) {
6d8df534 775 if(inindices[k]==id[l]) copyindex=kFALSE;
fee9ef0d 776 }
777 if(copyindex) {
778 outindices[j] = inindices[k]; j++;
779 }
780 }
4d023a8f 781 outVtx->SetIndices(nIndices,outindices);
fee9ef0d 782 delete [] outindices;
783
919e537f 784 /*
785 printf("Vertex before removing tracks:");
fee9ef0d 786 inVtx->PrintStatus();
787 inVtx->PrintIndices();
919e537f 788 printf("Vertex after removing tracks:");
fee9ef0d 789 outVtx->PrintStatus();
790 outVtx->PrintIndices();
919e537f 791 */
fee9ef0d 792
793 return outVtx;
794}
795//---------------------------------------------------------------------------
a00021a7 796void AliVertexerTracks::SetCuts(Double_t *cuts)
797{
798//
799// Cut values
800//
801 SetDCAcut(cuts[0]);
802 SetDCAcutIter0(cuts[1]);
803 SetMaxd0z0(cuts[2]);
5b78b791 804 if(fMode==0 && cuts[3]<0) SetITSrefitNotRequired();
805 SetMinClusters((Int_t)(TMath::Abs(cuts[3])));
a00021a7 806 SetMinTracks((Int_t)(cuts[4]));
807 SetNSigmad0(cuts[5]);
808 SetMinDetFitter(cuts[6]);
809 SetMaxTgl(cuts[7]);
810 SetFiducialRZ(cuts[8],cuts[9]);
811
812 return;
813}
814//---------------------------------------------------------------------------
f09c879d 815void AliVertexerTracks::SetITSMode(Double_t dcacut,
816 Double_t dcacutIter0,
817 Double_t maxd0z0,
6b29d399 818 Int_t minCls,
f09c879d 819 Int_t mintrks,
820 Double_t nsigma,
821 Double_t mindetfitter,
6b29d399 822 Double_t maxtgl,
823 Double_t fidR,
824 Double_t fidZ)
f09c879d 825{
826//
827// Cut values for ITS mode
828//
829 fMode = 0;
5b78b791 830 if(minCls>0) {
831 SetITSrefitRequired();
832 } else {
833 SetITSrefitNotRequired();
834 }
f09c879d 835 SetDCAcut(dcacut);
836 SetDCAcutIter0(dcacutIter0);
837 SetMaxd0z0(maxd0z0);
5b78b791 838 SetMinClusters(TMath::Abs(minCls));
f09c879d 839 SetMinTracks(mintrks);
840 SetNSigmad0(nsigma);
841 SetMinDetFitter(mindetfitter);
842 SetMaxTgl(maxtgl);
6b29d399 843 SetFiducialRZ(fidR,fidZ);
f09c879d 844
845 return;
846}
847//---------------------------------------------------------------------------
848void AliVertexerTracks::SetTPCMode(Double_t dcacut,
849 Double_t dcacutIter0,
850 Double_t maxd0z0,
6b29d399 851 Int_t minCls,
f09c879d 852 Int_t mintrks,
853 Double_t nsigma,
854 Double_t mindetfitter,
6b29d399 855 Double_t maxtgl,
856 Double_t fidR,
857 Double_t fidZ)
f09c879d 858{
859//
860// Cut values for TPC mode
861//
862 fMode = 1;
863 SetITSrefitNotRequired();
864 SetDCAcut(dcacut);
865 SetDCAcutIter0(dcacutIter0);
866 SetMaxd0z0(maxd0z0);
6b29d399 867 SetMinClusters(minCls);
f09c879d 868 SetMinTracks(mintrks);
869 SetNSigmad0(nsigma);
870 SetMinDetFitter(mindetfitter);
871 SetMaxTgl(maxtgl);
6b29d399 872 SetFiducialRZ(fidR,fidZ);
f09c879d 873
874 return;
875}
876//---------------------------------------------------------------------------
fee9ef0d 877void AliVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped)
878{
879//
880// Mark the tracks not to be used in the vertex reconstruction.
881// Tracks are identified by AliESDtrack::GetID()
882//
883 fNTrksToSkip = n; fTrksToSkip = new Int_t[n];
50ff0bcd 884 for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
885 return;
886}
887//---------------------------------------------------------------------------
fee9ef0d 888void AliVertexerTracks::SetVtxStart(AliESDVertex *vtx)
889{
50ff0bcd 890//
891// Set initial vertex knowledge
892//
893 vtx->GetXYZ(fNominalPos);
894 vtx->GetCovMatrix(fNominalCov);
fee9ef0d 895 SetConstraintOn();
50ff0bcd 896 return;
897}
2d57349e 898//---------------------------------------------------------------------------
fee9ef0d 899void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
dc3719f3 900{
6d8df534 901 AliExternalTrackParam *track1;
6d8df534 902 const Int_t knacc = (Int_t)fTrkArraySel.GetEntriesFast();
f54bad66 903 static TClonesArray linarray("AliStrLine",knacc);
dc3719f3 904 for(Int_t i=0; i<knacc; i++){
6d8df534 905 track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
dc3719f3 906 Double_t alpha=track1->GetAlpha();
907 Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
908 Double_t pos[3],dir[3],sigmasq[3];
f09c879d 909 track1->GetXYZAt(mindist,GetFieldkG(),pos);
910 track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
dc3719f3 911 sigmasq[0]=TMath::Sin(alpha)*TMath::Sin(alpha)*track1->GetSigmaY2();
912 sigmasq[1]=TMath::Cos(alpha)*TMath::Cos(alpha)*track1->GetSigmaY2();
913 sigmasq[2]=track1->GetSigmaZ2();
146c29df 914 TMatrixD ri(3,1);
915 TMatrixD wWi(3,3);
916 if(!TrackToPoint(track1,ri,wWi)) continue;
917 Double_t wmat[9];
918 Int_t iel=0;
919 for(Int_t ia=0;ia<3;ia++){
920 for(Int_t ib=0;ib<3;ib++){
921 wmat[iel]=wWi(ia,ib);
922 iel++;
923 }
924 }
f54bad66 925 new(linarray[i]) AliStrLine(pos,sigmasq,wmat,dir);
dc3719f3 926 }
f54bad66 927 fVert=TrackletVertexFinder(&linarray,optUseWeights);
928 linarray.Clear("C");
dc3719f3 929}
930//---------------------------------------------------------------------------
146c29df 931AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t optUseWeights)
fee9ef0d 932{
2d57349e 933 // Calculate the point at minimum distance to prepared tracks
934
dc3719f3 935 const Int_t knacc = (Int_t)lines->GetEntriesFast();
936 Double_t initPos[3]={0.,0.,0.};
2d57349e 937
2d57349e 938 Double_t (*vectP0)[3]=new Double_t [knacc][3];
939 Double_t (*vectP1)[3]=new Double_t [knacc][3];
940
941 Double_t sum[3][3];
942 Double_t dsum[3]={0,0,0};
146c29df 943 TMatrixD sumWi(3,3);
944 for(Int_t i=0;i<3;i++){
945 for(Int_t j=0;j<3;j++){
946 sum[i][j]=0;
947 sumWi(i,j)=0.;
948 }
949 }
950
2d57349e 951 for(Int_t i=0; i<knacc; i++){
f54bad66 952 AliStrLine* line1 = (AliStrLine*)lines->At(i);
dc3719f3 953 Double_t p0[3],cd[3],sigmasq[3];
146c29df 954 Double_t wmat[9];
2d57349e 955 line1->GetP0(p0);
956 line1->GetCd(cd);
dc3719f3 957 line1->GetSigma2P0(sigmasq);
146c29df 958 line1->GetWMatrix(wmat);
959 TMatrixD wWi(3,3);
960 Int_t iel=0;
961 for(Int_t ia=0;ia<3;ia++){
962 for(Int_t ib=0;ib<3;ib++){
963 wWi(ia,ib)=wmat[iel];
964 iel++;
965 }
966 }
967
968 sumWi+=wWi;
969
2d57349e 970 Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
971 vectP0[i][0]=p0[0];
972 vectP0[i][1]=p0[1];
973 vectP0[i][2]=p0[2];
974 vectP1[i][0]=p1[0];
975 vectP1[i][1]=p1[1];
976 vectP1[i][2]=p1[2];
977
978 Double_t matr[3][3];
979 Double_t dknow[3];
980 if(optUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
dc3719f3 981 else GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
982
2d57349e 983
984 for(Int_t iii=0;iii<3;iii++){
985 dsum[iii]+=dknow[iii];
986 for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
987 }
2d57349e 988 }
146c29df 989
990 TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
991 Double_t covmatrix[6];
992 covmatrix[0] = invsumWi(0,0);
993 covmatrix[1] = invsumWi(0,1);
994 covmatrix[2] = invsumWi(1,1);
995 covmatrix[3] = invsumWi(0,2);
996 covmatrix[4] = invsumWi(1,2);
997 covmatrix[5] = invsumWi(2,2);
998
2d57349e 999 Double_t vett[3][3];
1000 Double_t det=GetDeterminant3X3(sum);
dc3719f3 1001 Double_t sigma=0;
2d57349e 1002
fee9ef0d 1003 if(det!=0){
1004 for(Int_t zz=0;zz<3;zz++){
1005 for(Int_t ww=0;ww<3;ww++){
1006 for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
1007 }
1008 for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
1009 initPos[zz]=GetDeterminant3X3(vett)/det;
1010 }
1011
1012
1013 for(Int_t i=0; i<knacc; i++){
1014 Double_t p0[3]={0,0,0},p1[3]={0,0,0};
1015 for(Int_t ii=0;ii<3;ii++){
1016 p0[ii]=vectP0[i][ii];
1017 p1[ii]=vectP1[i][ii];
1018 }
1019 sigma+=GetStrLinMinDist(p0,p1,initPos);
1020 }
1021
f09c879d 1022 if(sigma>0.) {sigma=TMath::Sqrt(sigma);}else{sigma=999;}
fee9ef0d 1023 }else{
50ff0bcd 1024 sigma=999;
1025 }
146c29df 1026 AliESDVertex theVert(initPos,covmatrix,99999.,knacc);
1027 theVert.SetDispersion(sigma);
50ff0bcd 1028 delete vectP0;
1029 delete vectP1;
dc3719f3 1030 return theVert;
50ff0bcd 1031}
1032//---------------------------------------------------------------------------
6d8df534 1033Bool_t AliVertexerTracks::TrackToPoint(AliExternalTrackParam *t,
0bf2e2b4 1034 TMatrixD &ri,TMatrixD &wWi,
1035 Bool_t uUi3by3) const
fee9ef0d 1036{
1037//
6d8df534 1038// Extract from the AliExternalTrackParam the global coordinates ri and covariance matrix
fee9ef0d 1039// wWi of the space point that it represents (to be used in VertexFitter())
1040//
1041
1042
1043 Double_t rotAngle = t->GetAlpha();
1044 if(rotAngle<0.) rotAngle += 2.*TMath::Pi();
1045 Double_t cosRot = TMath::Cos(rotAngle);
1046 Double_t sinRot = TMath::Sin(rotAngle);
1047
1048 ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
1049 ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
1050 ri(2,0) = t->GetZ();
1051
0bf2e2b4 1052
1053
1054 if(!uUi3by3) {
1055 // matrix to go from global (x,y,z) to local (y,z);
1056 TMatrixD qQi(2,3);
1057 qQi(0,0) = -sinRot;
1058 qQi(0,1) = cosRot;
1059 qQi(0,2) = 0.;
1060 qQi(1,0) = 0.;
1061 qQi(1,1) = 0.;
1062 qQi(1,2) = 1.;
1063
1064 // covariance matrix of local (y,z) - inverted
0bf2e2b4 1065 TMatrixD uUi(2,2);
6d8df534 1066 uUi(0,0) = t->GetSigmaY2();
1067 uUi(0,1) = t->GetSigmaZY();
1068 uUi(1,0) = t->GetSigmaZY();
1069 uUi(1,1) = t->GetSigmaZ2();
919e537f 1070 //printf(" Ui :");
1071 //printf(" %f %f",uUi(0,0),uUi(0,1));
1072 //printf(" %f %f",uUi(1,0),uUi(1,1));
0bf2e2b4 1073
1074 if(uUi.Determinant() <= 0.) return kFALSE;
1075 TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1076
1077 // weights matrix: wWi = qQiT * uUiInv * qQi
1078 TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1079 TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1080 wWi = m;
1081 } else {
1082 if(fVert.GetNContributors()<1) AliFatal("Vertex from finder is empty");
1083 // matrix to go from global (x,y,z) to local (x,y,z);
1084 TMatrixD qQi(3,3);
1085 qQi(0,0) = cosRot;
1086 qQi(0,1) = sinRot;
1087 qQi(0,2) = 0.;
1088 qQi(1,0) = -sinRot;
1089 qQi(1,1) = cosRot;
1090 qQi(1,2) = 0.;
1091 qQi(2,0) = 0.;
1092 qQi(2,1) = 0.;
1093 qQi(2,2) = 1.;
1094
1095 // covariance of fVert along the track
1096 Double_t p[3],pt,ptot;
1097 t->GetPxPyPz(p);
1098 pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
1099 ptot = TMath::Sqrt(pt*pt+p[2]*p[2]);
3b768a1d 1100 Double_t cphi = p[0]/pt; //cos(phi)=px/pt
1101 Double_t sphi = p[1]/pt; //sin(phi)=py/pt
1102 Double_t clambda = pt/ptot; //cos(lambda)=pt/ptot
1103 Double_t slambda = p[2]/ptot; //sin(lambda)=pz/ptot
0bf2e2b4 1104 Double_t covfVert[6];
1105 fVert.GetCovMatrix(covfVert);
1106 Double_t covfVertalongt =
3b768a1d 1107 covfVert[0]*cphi*cphi*clambda*clambda
1108 +covfVert[1]*2.*cphi*sphi*clambda*clambda
1109 +covfVert[3]*2.*cphi*clambda*slambda
1110 +covfVert[2]*sphi*sphi*clambda*clambda
1111 +covfVert[4]*2.*sphi*clambda*slambda
1112 +covfVert[5]*slambda*slambda;
0bf2e2b4 1113 // covariance matrix of local (x,y,z) - inverted
1114 TMatrixD uUi(3,3);
3b768a1d 1115 uUi(0,0) = covfVertalongt * fnSigmaForUi00 * fnSigmaForUi00;
919e537f 1116 AliDebug(1,Form("=====> sqrtUi00 cm %f",TMath::Sqrt(uUi(0,0))));
0bf2e2b4 1117 uUi(0,1) = 0.;
1118 uUi(0,2) = 0.;
1119 uUi(1,0) = 0.;
6d8df534 1120 uUi(1,1) = t->GetSigmaY2();
1121 uUi(1,2) = t->GetSigmaZY();
0bf2e2b4 1122 uUi(2,0) = 0.;
6d8df534 1123 uUi(2,1) = t->GetSigmaZY();
1124 uUi(2,2) = t->GetSigmaZ2();
0bf2e2b4 1125 //printf(" Ui :\n");
1126 //printf(" %f %f\n",uUi(0,0),uUi(0,1));
1127 //printf(" %f %f\n",uUi(1,0),uUi(1,1));
fee9ef0d 1128
0bf2e2b4 1129 if(uUi.Determinant() <= 0.) return kFALSE;
1130 TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1131
1132 // weights matrix: wWi = qQiT * uUiInv * qQi
1133 TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1134 TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1135 wWi = m;
1136 }
1137
fee9ef0d 1138
1139 return kTRUE;
1140}
1141//---------------------------------------------------------------------------
6d8df534 1142void AliVertexerTracks::TooFewTracks()
fee9ef0d 1143{
50ff0bcd 1144//
6d8df534 1145// When the number of tracks is < fMinTracks,
1146// deal with vertices not found and prepare to exit
50ff0bcd 1147//
919e537f 1148 AliDebug(1,"TooFewTracks");
2d57349e 1149
50ff0bcd 1150 Double_t pos[3],err[3];
50ff0bcd 1151 pos[0] = fNominalPos[0];
1152 err[0] = TMath::Sqrt(fNominalCov[0]);
1153 pos[1] = fNominalPos[1];
1154 err[1] = TMath::Sqrt(fNominalCov[2]);
6d8df534 1155 pos[2] = fNominalPos[2];
1156 err[2] = TMath::Sqrt(fNominalCov[5]);
1157 Int_t ncontr = (err[0]>1. ? -1 : -3);
5b78b791 1158 if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
50ff0bcd 1159 fCurrentVertex = new AliESDVertex(pos,err);
1160 fCurrentVertex->SetNContributors(ncontr);
2d57349e 1161
fee9ef0d 1162 if(fConstraint) {
1163 fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
1164 } else {
50ff0bcd 1165 fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
fee9ef0d 1166 }
2d57349e 1167
6d8df534 1168 if(!fTrkArraySel.IsEmpty()) fTrkArraySel.Delete();
1169 if(fIdSel) {delete [] fIdSel; fIdSel=NULL;}
1170 if(fTrksToSkip) {delete [] fTrksToSkip; fTrksToSkip=NULL;}
1171
50ff0bcd 1172 return;
1173}
1174//---------------------------------------------------------------------------
fee9ef0d 1175void AliVertexerTracks::VertexFinder(Int_t optUseWeights)
1176{
2d57349e 1177
50ff0bcd 1178 // Get estimate of vertex position in (x,y) from tracks DCA
1179
1180 Double_t initPos[3];
1181 initPos[2] = 0.;
1182 for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
6d8df534 1183 Int_t nacc = (Int_t)fTrkArraySel.GetEntriesFast();
50ff0bcd 1184 Double_t aver[3]={0.,0.,0.};
1185 Double_t aversq[3]={0.,0.,0.};
1186 Double_t sigmasq[3]={0.,0.,0.};
1187 Double_t sigma=0;
1188 Int_t ncombi = 0;
6d8df534 1189 AliExternalTrackParam *track1;
1190 AliExternalTrackParam *track2;
50ff0bcd 1191 Double_t pos[3],dir[3];
1192 Double_t alpha,mindist;
2d57349e 1193
50ff0bcd 1194 for(Int_t i=0; i<nacc; i++){
6d8df534 1195 track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
50ff0bcd 1196 alpha=track1->GetAlpha();
1197 mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
f09c879d 1198 track1->GetXYZAt(mindist,GetFieldkG(),pos);
1199 track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
50ff0bcd 1200 AliStrLine *line1 = new AliStrLine(pos,dir);
2d57349e 1201
50ff0bcd 1202 // AliStrLine *line1 = new AliStrLine();
f09c879d 1203 // track1->ApproximateHelixWithLine(mindist,GetFieldkG(),line1);
50ff0bcd 1204
1205 for(Int_t j=i+1; j<nacc; j++){
6d8df534 1206 track2 = (AliExternalTrackParam*)fTrkArraySel.At(j);
50ff0bcd 1207 alpha=track2->GetAlpha();
1208 mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
f09c879d 1209 track2->GetXYZAt(mindist,GetFieldkG(),pos);
1210 track2->GetPxPyPzAt(mindist,GetFieldkG(),dir);
50ff0bcd 1211 AliStrLine *line2 = new AliStrLine(pos,dir);
1212 // AliStrLine *line2 = new AliStrLine();
f09c879d 1213 // track2->ApproximateHelixWithLine(mindist,GetFieldkG(),line2);
50ff0bcd 1214 Double_t distCA=line2->GetDCA(line1);
fee9ef0d 1215 //printf("%d %d %f\n",i,j,distCA);
1216 if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
50ff0bcd 1217 Double_t pnt1[3],pnt2[3],crosspoint[3];
2d57349e 1218
50ff0bcd 1219 if(optUseWeights<=0){
1220 Int_t retcode = line2->Cross(line1,crosspoint);
1221 if(retcode>=0){
1222 ncombi++;
1223 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1224 for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1225 }
1226 }
1227 if(optUseWeights>0){
1228 Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2);
1229 if(retcode>=0){
5c03260f 1230 Double_t cs, sn;
50ff0bcd 1231 alpha=track1->GetAlpha();
1232 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1233 Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
1234 alpha=track2->GetAlpha();
1235 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1236 Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
1237 Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
1238 Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
1239 Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
1240 Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
1241 crosspoint[0]=wx1*pnt1[0] + wx2*pnt2[0];
1242 crosspoint[1]=wy1*pnt1[1] + wy2*pnt2[1];
1243 crosspoint[2]=wz1*pnt1[2] + wz2*pnt2[2];
1244
1245 ncombi++;
1246 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1247 for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1248 }
1249 }
1250 }
1251 delete line2;
1252 }
1253 delete line1;
71d84967 1254 }
50ff0bcd 1255 if(ncombi>0){
1256 for(Int_t jj=0;jj<3;jj++){
1257 initPos[jj] = aver[jj]/ncombi;
fee9ef0d 1258 //printf("%f\n",initPos[jj]);
50ff0bcd 1259 aversq[jj]/=ncombi;
1260 sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
1261 sigma+=sigmasq[jj];
1262 }
1263 sigma=TMath::Sqrt(TMath::Abs(sigma));
1264 }
1265 else {
1266 Warning("VertexFinder","Finder did not succed");
1267 sigma=999;
1268 }
1269 fVert.SetXYZ(initPos);
1270 fVert.SetDispersion(sigma);
1271 fVert.SetNContributors(ncombi);
2d57349e 1272}
ec1be5d5 1273//---------------------------------------------------------------------------
6d8df534 1274void AliVertexerTracks::VertexFitter()
fee9ef0d 1275{
ec1be5d5 1276//
1277// The optimal estimate of the vertex position is given by a "weighted
0bf2e2b4 1278// average of tracks positions".
fee9ef0d 1279// Original method: V. Karimaki, CMS Note 97/0051
ec1be5d5 1280//
6d8df534 1281 Bool_t useConstraint = fConstraint;
ec1be5d5 1282 Double_t initPos[3];
0bf2e2b4 1283 if(!fOnlyFitter) {
1284 fVert.GetXYZ(initPos);
1285 } else {
1286 initPos[0]=fNominalPos[0];
1287 initPos[1]=fNominalPos[1];
1288 initPos[2]=fNominalPos[2];
1289 }
1290
6d8df534 1291 Int_t nTrksSel = (Int_t)fTrkArraySel.GetEntries();
1292 if(nTrksSel==1) useConstraint=kTRUE;
919e537f 1293 AliDebug(1,"--- VertexFitter(): start");
1294 AliDebug(1,Form(" Number of tracks in array: %d\n",nTrksSel));
1295 AliDebug(1,Form(" Minimum # tracks required in fit: %d\n",fMinTracks));
1296 AliDebug(1,Form(" Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]));
1297 if(useConstraint) AliDebug(1,Form(" This vertex will be used in fit: (%f+-%f,%f+-%f)\n",fNominalPos[0],TMath::Sqrt(fNominalCov[0]),fNominalPos[1],TMath::Sqrt(fNominalCov[2])));
ec1be5d5 1298
0bf2e2b4 1299 // special treatment for few-tracks fits (e.g. secondary vertices)
6d8df534 1300 Bool_t uUi3by3 = kFALSE; if(nTrksSel<5 && !useConstraint) uUi3by3 = kTRUE;
0bf2e2b4 1301
ec1be5d5 1302 Int_t i,j,k,step=0;
1303 TMatrixD rv(3,1);
1304 TMatrixD vV(3,3);
1305 rv(0,0) = initPos[0];
1306 rv(1,0) = initPos[1];
1307 rv(2,0) = 0.;
1308 Double_t xlStart,alpha;
6d8df534 1309 Int_t nTrksUsed;
fee9ef0d 1310 Double_t chi2,chi2i,chi2b;
6d8df534 1311 AliExternalTrackParam *t = 0;
ec1be5d5 1312 Int_t failed = 0;
1313
50ff0bcd 1314 // initial vertex covariance matrix
1315 TMatrixD vVb(3,3);
1316 vVb(0,0) = fNominalCov[0];
1317 vVb(0,1) = fNominalCov[1];
1318 vVb(0,2) = 0.;
1319 vVb(1,0) = fNominalCov[1];
1320 vVb(1,1) = fNominalCov[2];
1321 vVb(1,2) = 0.;
1322 vVb(2,0) = 0.;
1323 vVb(2,1) = 0.;
1324 vVb(2,2) = fNominalCov[5];
1325 TMatrixD vVbInv(TMatrixD::kInverted,vVb);
1326 TMatrixD rb(3,1);
1327 rb(0,0) = fNominalPos[0];
1328 rb(1,0) = fNominalPos[1];
1329 rb(2,0) = fNominalPos[2];
1330 TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb);
1331
1332
eae1677e 1333 // 2 steps:
1334 // 1st - estimate of vtx using all tracks
1335 // 2nd - estimate of global chi2
1336 for(step=0; step<2; step++) {
919e537f 1337 AliDebug(1,Form(" step = %d\n",step));
ec1be5d5 1338 chi2 = 0.;
6d8df534 1339 nTrksUsed = 0;
ec1be5d5 1340
0bf2e2b4 1341 if(step==1) { initPos[0]=rv(0,0); initPos[0]=rv(1,0); }
1342
ec1be5d5 1343 TMatrixD sumWiri(3,1);
1344 TMatrixD sumWi(3,3);
1345 for(i=0; i<3; i++) {
1346 sumWiri(i,0) = 0.;
1347 for(j=0; j<3; j++) sumWi(j,i) = 0.;
1348 }
1349
fee9ef0d 1350 // mean vertex constraint
1351 if(useConstraint) {
50ff0bcd 1352 for(i=0;i<3;i++) {
1353 sumWiri(i,0) += vVbInvrb(i,0);
1354 for(k=0;k<3;k++) sumWi(i,k) += vVbInv(i,k);
07680cae 1355 }
fee9ef0d 1356 // chi2
1357 TMatrixD deltar = rv; deltar -= rb;
1358 TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar);
1359 chi2b = deltar(0,0)*vVbInvdeltar(0,0)+
1360 deltar(1,0)*vVbInvdeltar(1,0)+
1361 deltar(2,0)*vVbInvdeltar(2,0);
1362 chi2 += chi2b;
07680cae 1363 }
1364
ec1be5d5 1365 // loop on tracks
6d8df534 1366 for(k=0; k<nTrksSel; k++) {
ec1be5d5 1367 // get track from track array
6d8df534 1368 t = (AliExternalTrackParam*)fTrkArraySel.At(k);
ec1be5d5 1369 alpha = t->GetAlpha();
1370 xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
fee9ef0d 1371 // to vtxSeed (from finder)
6d8df534 1372 t->PropagateTo(xlStart,GetFieldkG());
fee9ef0d 1373
ec1be5d5 1374 // vector of track global coordinates
1375 TMatrixD ri(3,1);
fee9ef0d 1376 // covariance matrix of ri
1377 TMatrixD wWi(3,3);
1378
1379 // get space point from track
0bf2e2b4 1380 if(!TrackToPoint(t,ri,wWi,uUi3by3)) continue;
fee9ef0d 1381 TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
ec1be5d5 1382
1383 // track chi2
1384 TMatrixD deltar = rv; deltar -= ri;
1385 TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
1386 chi2i = deltar(0,0)*wWideltar(0,0)+
1387 deltar(1,0)*wWideltar(1,0)+
1388 deltar(2,0)*wWideltar(2,0);
1389
ec1be5d5 1390 // add to total chi2
1391 chi2 += chi2i;
1392
ec1be5d5 1393 sumWiri += wWiri;
1394 sumWi += wWi;
1395
6d8df534 1396 nTrksUsed++;
ec1be5d5 1397 } // end loop on tracks
1398
6d8df534 1399 if(nTrksUsed < fMinTracks) {
ec1be5d5 1400 failed=1;
1401 continue;
1402 }
1403
1404 Double_t determinant = sumWi.Determinant();
f09c879d 1405 if(determinant < fMinDetFitter) {
919e537f 1406 AliDebug(1,Form("det(V) = %f (<%f)\n",determinant,fMinDetFitter));
ec1be5d5 1407 failed=1;
1408 continue;
1409 }
1410
0bf2e2b4 1411 if(step==0) {
1412 // inverted of weights matrix
1413 TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1414 vV = invsumWi;
1415 // position of primary vertex
1416 rv.Mult(vV,sumWiri);
1417 }
eae1677e 1418 } // end loop on the 2 steps
ec1be5d5 1419
ec1be5d5 1420 if(failed) {
6d8df534 1421 TooFewTracks();
ec1be5d5 1422 return;
1423 }
1424
1425 Double_t position[3];
1426 position[0] = rv(0,0);
1427 position[1] = rv(1,0);
1428 position[2] = rv(2,0);
1429 Double_t covmatrix[6];
1430 covmatrix[0] = vV(0,0);
1431 covmatrix[1] = vV(0,1);
1432 covmatrix[2] = vV(1,1);
1433 covmatrix[3] = vV(0,2);
1434 covmatrix[4] = vV(1,2);
1435 covmatrix[5] = vV(2,2);
1436
dc3719f3 1437 // for correct chi2/ndf, count constraint as additional "track"
6d8df534 1438 if(fConstraint) nTrksUsed++;
ec1be5d5 1439 // store data in the vertex object
5b78b791 1440 if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
6d8df534 1441 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
ec1be5d5 1442
919e537f 1443 AliDebug(1," Vertex after fit:");
1444 AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
1445 AliDebug(1,"--- VertexFitter(): finish\n");
1446
ec1be5d5 1447
1448 return;
1449}
50ff0bcd 1450//----------------------------------------------------------------------------
6d8df534 1451AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray,
1452 UShort_t *id,
1453 Bool_t optUseFitter,
1454 Bool_t optPropagate)
fee9ef0d 1455{
eae1677e 1456//
6d8df534 1457// Return vertex from tracks (AliExternalTrackParam) in array
eae1677e 1458//
5b78b791 1459 fCurrentVertex = 0;
fee9ef0d 1460 SetConstraintOff();
1461
50ff0bcd 1462 // get tracks and propagate them to initial vertex position
6d8df534 1463 fIdSel = new UShort_t[(Int_t)trkArray->GetEntriesFast()];
1464 Int_t nTrksSel = PrepareTracks(*trkArray,id,0);
1465 if(nTrksSel < TMath::Max(2,fMinTracks)) {
1466 TooFewTracks();
fee9ef0d 1467 return fCurrentVertex;
50ff0bcd 1468 }
1469
6d8df534 1470 // vertex finder
fee9ef0d 1471 switch (fAlgo) {
1472 case 1: StrLinVertexFinderMinDist(1); break;
1473 case 2: StrLinVertexFinderMinDist(0); break;
1474 case 3: HelixVertexFinder(); break;
1475 case 4: VertexFinder(1); break;
1476 case 5: VertexFinder(0); break;
1477 default: printf("Wrong algorithm\n"); break;
1478 }
919e537f 1479 AliDebug(1," Vertex finding completed\n");
fee9ef0d 1480
1481 // vertex fitter
6d8df534 1482 if(optUseFitter) {
1483 VertexFitter();
1484 } else {
fee9ef0d 1485 Double_t position[3]={fVert.GetXv(),fVert.GetYv(),fVert.GetZv()};
146c29df 1486 Double_t covmatrix[6];
1487 fVert.GetCovMatrix(covmatrix);
fee9ef0d 1488 Double_t chi2=99999.;
6d8df534 1489 Int_t nTrksUsed=fVert.GetNContributors();
1490 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
fee9ef0d 1491 }
1492 fCurrentVertex->SetDispersion(fVert.GetDispersion());
1493
1494
1495 // set indices of used tracks and propagate track to found vertex
50ff0bcd 1496 UShort_t *indices = 0;
6d8df534 1497 Double_t d0z0[2],covd0z0[3];
1498 AliExternalTrackParam *t = 0;
fee9ef0d 1499 if(fCurrentVertex->GetNContributors()>0) {
1500 indices = new UShort_t[fCurrentVertex->GetNContributors()];
6d8df534 1501 for(Int_t jj=0; jj<(Int_t)fTrkArraySel.GetEntriesFast(); jj++) {
1502 indices[jj] = fIdSel[jj];
1503 t = (AliExternalTrackParam*)fTrkArraySel.At(jj);
1504 if(optPropagate && optUseFitter) {
fee9ef0d 1505 if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
6d8df534 1506 t->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
919e537f 1507 AliDebug(1,Form("Track %d propagated to found vertex",jj));
6d8df534 1508 } else {
fee9ef0d 1509 AliWarning("Found vertex outside beam pipe!");
1510 }
1511 }
50ff0bcd 1512 }
fee9ef0d 1513 fCurrentVertex->SetIndices(fCurrentVertex->GetNContributors(),indices);
50ff0bcd 1514 }
6d8df534 1515
1516 // clean up
1517 delete [] indices; indices=NULL;
1518 delete [] fIdSel; fIdSel=NULL;
1519 fTrkArraySel.Delete();
50ff0bcd 1520
fee9ef0d 1521 return fCurrentVertex;
eae1677e 1522}
50ff0bcd 1523//----------------------------------------------------------------------------
6d8df534 1524AliESDVertex* AliVertexerTracks::VertexForSelectedESDTracks(TObjArray *trkArray,
1525 Bool_t optUseFitter,
1526 Bool_t optPropagate)
fee9ef0d 1527{
50ff0bcd 1528//
6d8df534 1529// Return vertex from array of ESD tracks
50ff0bcd 1530//
5834e9a1 1531
6d8df534 1532 Int_t nTrks = (Int_t)trkArray->GetEntriesFast();
1533 UShort_t *id = new UShort_t[nTrks];
1534
1535 AliESDtrack *esdt = 0;
50ff0bcd 1536 for(Int_t i=0; i<nTrks; i++){
6d8df534 1537 esdt = (AliESDtrack*)trkArray->At(i);
1538 id[i] = (UShort_t)esdt->GetID();
50ff0bcd 1539 }
1540
6d8df534 1541 VertexForSelectedTracks(trkArray,id,optUseFitter,optPropagate);
1542
1543 delete [] id; id=NULL;
1544
1545 return fCurrentVertex;
50ff0bcd 1546}
1547//--------------------------------------------------------------------------