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