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