]> 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
Adding TString in the included files.
[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
d993e038 603 Bool_t propagateOK = kFALSE, cutond0z0 = kTRUE;
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]);
d993e038 611 AliDebug(1,Form("normdistx %f %f %f",fCurrentVertex->GetXv(),fNominalPos[0],TMath::Sqrt(sigmaCurr[0]*sigmaCurr[0]+fNominalCov[0])));
612 AliDebug(1,Form("normdisty %f %f %f",fCurrentVertex->GetYv(),fNominalPos[1],TMath::Sqrt(sigmaCurr[1]*sigmaCurr[1]+fNominalCov[2])));
613 AliDebug(1,Form("sigmaCurr %f %f %f",sigmaCurr[0],sigmaCurr[1],TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2])));
fee9ef0d 614 if(normdistx < 3. && normdisty < 3. &&
80d36431 615 (sigmaCurr[0]+sigmaCurr[1])<(TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2]))) {
18eba197 616 propagateOK = track->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
50ff0bcd 617 } else {
18eba197 618 propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
d993e038 619 if(fConstraint) cutond0z0=kFALSE;
50ff0bcd 620 }
621 }
622
18eba197 623 if(!propagateOK) {
919e537f 624 AliDebug(1," rejected");
18eba197 625 delete track; continue;
626 }
627
43c9dae1 628 sigmad0 = TMath::Sqrt(covd0z0[0]);
629 maxd0rphi = fNSigma*sigmad0;
fee9ef0d 630 if(optImpParCut==1) maxd0rphi *= 5.;
6b29d399 631 maxd0rphi = TMath::Min(maxd0rphi,fFiducialR);
d993e038 632 //sigmad0z0 = TMath::Sqrt(covd0z0[0]+covd0z0[2]);
fee9ef0d 633
919e537f 634 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 635
636
637 //---- track selection based on impact parameters ----//
638
639 // always reject tracks outside fiducial volume
6b29d399 640 if(TMath::Abs(d0z0[0])>fFiducialR || TMath::Abs(d0z0[1])>fFiducialZ) {
919e537f 641 AliDebug(1," rejected");
f09c879d 642 delete track; continue;
643 }
644
645 // during iterations 1 and 2 , reject tracks with d0rphi > maxd0rphi
646 if(optImpParCut>0 && TMath::Abs(d0z0[0])>maxd0rphi) {
919e537f 647 AliDebug(1," rejected");
f09c879d 648 delete track; continue;
649 }
fee9ef0d 650
43c9dae1 651 // if fConstraint=kFALSE, during iterations 1 and 2 ||
652 // if fConstraint=kTRUE, during iteration 2,
f09c879d 653 // select tracks with d0oplusz0 < fMaxd0z0
43c9dae1 654 if((!fConstraint && optImpParCut>0 && fVert.GetNContributors()>0) ||
d993e038 655 ( fConstraint && optImpParCut==2 && cutond0z0)) {
43c9dae1 656 if(nTrksOrig>=3 &&
f09c879d 657 TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1])>fMaxd0z0) {
919e537f 658 AliDebug(1," rejected");
fee9ef0d 659 delete track; continue;
660 }
661 }
43c9dae1 662
f09c879d 663 // track passed all selections
6d8df534 664 fTrkArraySel.AddLast(track);
665 fIdSel[nTrksSel] = idOrig[i];
666 nTrksSel++;
f09c879d 667 } // end loop on tracks
50ff0bcd 668
669 delete initVertex;
670
6d8df534 671 return nTrksSel;
50ff0bcd 672}
f09c879d 673//----------------------------------------------------------------------------
674Bool_t AliVertexerTracks::PropagateTrackTo(AliExternalTrackParam *track,
675 Double_t xToGo) {
676 //----------------------------------------------------------------
677 // COPIED from AliTracker
678 //
679 // Propagates the track to the plane X=xk (cm).
680 //
681 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
682 //----------------------------------------------------------------
683
684 const Double_t kEpsilon = 0.00001;
685 Double_t xpos = track->GetX();
686 Double_t dir = (xpos<xToGo) ? 1. : -1.;
687 Double_t maxStep = 5;
688 Double_t maxSnp = 0.8;
689 //
690 while ( (xToGo-xpos)*dir > kEpsilon){
691 Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
692 Double_t x = xpos+step;
693 Double_t xyz0[3],xyz1[3];
694 track->GetXYZ(xyz0); //starting global position
695
696 if(!track->GetXYZAt(x,GetFieldkG(),xyz1)) return kFALSE; // no prolongation
697 xyz1[2]+=kEpsilon; // waiting for bug correction in geo
698
699 if(TMath::Abs(track->GetSnpAt(x,GetFieldkG())) >= maxSnp) return kFALSE;
700 if(!track->PropagateTo(x,GetFieldkG())) return kFALSE;
701
702 if(TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
703 track->GetXYZ(xyz0); // global position
704 Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]);
705 //
706 Double_t ca=TMath::Cos(alphan-track->GetAlpha()),
707 sa=TMath::Sin(alphan-track->GetAlpha());
60e55aee 708 Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
f09c879d 709 Double_t sinNew = sf*ca - cf*sa;
710 if(TMath::Abs(sinNew) >= maxSnp) return kFALSE;
711 if(!track->Rotate(alphan)) return kFALSE;
712
713 xpos = track->GetX();
714 }
715 return kTRUE;
716}
50ff0bcd 717//---------------------------------------------------------------------------
fee9ef0d 718AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
6d8df534 719 TObjArray *trkArray,
720 UShort_t *id,
6b908072 721 Float_t *diamondxy) const
fee9ef0d 722{
50ff0bcd 723//
fee9ef0d 724// Removes tracks in trksTree from fit of inVtx
50ff0bcd 725//
fee9ef0d 726
146c29df 727 if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraint")) {
919e537f 728 printf("ERROR: primary vertex has no constraint: cannot remove tracks");
146c29df 729 return 0x0;
730 }
fee9ef0d 731 if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
919e537f 732 printf("WARNING: result of tracks' removal will be only approximately correct");
fee9ef0d 733
734 TMatrixD rv(3,1);
735 rv(0,0) = inVtx->GetXv();
736 rv(1,0) = inVtx->GetYv();
737 rv(2,0) = inVtx->GetZv();
738 TMatrixD vV(3,3);
739 Double_t cov[6];
740 inVtx->GetCovMatrix(cov);
741 vV(0,0) = cov[0];
742 vV(0,1) = cov[1]; vV(1,0) = cov[1];
743 vV(1,1) = cov[2];
744 vV(0,2) = cov[3]; vV(2,0) = cov[3];
745 vV(1,2) = cov[4]; vV(2,1) = cov[4];
746 vV(2,2) = cov[5];
747
748 TMatrixD sumWi(TMatrixD::kInverted,vV);
749 TMatrixD sumWiri(sumWi,TMatrixD::kMult,rv);
750
dad616ce 751 Int_t nUsedTrks = inVtx->GetNIndices();
fee9ef0d 752 Double_t chi2 = inVtx->GetChi2();
753
6d8df534 754 AliExternalTrackParam *track = 0;
755 Int_t ntrks = (Int_t)trkArray->GetEntriesFast();
fee9ef0d 756 for(Int_t i=0;i<ntrks;i++) {
6d8df534 757 track = (AliExternalTrackParam*)trkArray->At(i);
758 if(!inVtx->UsesTrack(id[i])) {
919e537f 759 printf("track %d was not used in vertex fit",id[i]);
fee9ef0d 760 continue;
761 }
762 Double_t alpha = track->GetAlpha();
763 Double_t xl = diamondxy[0]*TMath::Cos(alpha)+diamondxy[1]*TMath::Sin(alpha);
6d8df534 764 track->PropagateTo(xl,GetFieldkG());
fee9ef0d 765 // vector of track global coordinates
766 TMatrixD ri(3,1);
767 // covariance matrix of ri
768 TMatrixD wWi(3,3);
769
770 // get space point from track
771 if(!TrackToPoint(track,ri,wWi)) continue;
772
773 TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
774
775 sumWi -= wWi;
776 sumWiri -= wWiri;
777
6b908072 778 // track contribution to chi2
fee9ef0d 779 TMatrixD deltar = rv; deltar -= ri;
780 TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
781 Double_t chi2i = deltar(0,0)*wWideltar(0,0)+
782 deltar(1,0)*wWideltar(1,0)+
783 deltar(2,0)*wWideltar(2,0);
784 // remove from total chi2
785 chi2 -= chi2i;
786
787 nUsedTrks--;
788 if(nUsedTrks<2) {
919e537f 789 printf("Trying to remove too many tracks!");
fee9ef0d 790 return 0x0;
791 }
792 }
793
794 TMatrixD rvnew(3,1);
795 TMatrixD vVnew(3,3);
796
797 // new inverted of weights matrix
798 TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
799 vVnew = invsumWi;
800 // new position of primary vertex
801 rvnew.Mult(vVnew,sumWiri);
802
803 Double_t position[3];
804 position[0] = rvnew(0,0);
805 position[1] = rvnew(1,0);
806 position[2] = rvnew(2,0);
807 cov[0] = vVnew(0,0);
808 cov[1] = vVnew(0,1);
809 cov[2] = vVnew(1,1);
810 cov[3] = vVnew(0,2);
811 cov[4] = vVnew(1,2);
812 cov[5] = vVnew(2,2);
813
814 // store data in the vertex object
6b908072 815 AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,nUsedTrks+1); // the +1 is for the constraint
fee9ef0d 816 outVtx->SetTitle(inVtx->GetTitle());
fee9ef0d 817 UShort_t *inindices = inVtx->GetIndices();
dad616ce 818 Int_t nIndices = nUsedTrks;
4d023a8f 819 UShort_t *outindices = new UShort_t[nIndices];
fee9ef0d 820 Int_t j=0;
fee9ef0d 821 for(Int_t k=0; k<inVtx->GetNIndices(); k++) {
6d8df534 822 Bool_t copyindex=kTRUE;
fee9ef0d 823 for(Int_t l=0; l<ntrks; l++) {
6d8df534 824 if(inindices[k]==id[l]) copyindex=kFALSE;
fee9ef0d 825 }
826 if(copyindex) {
827 outindices[j] = inindices[k]; j++;
828 }
829 }
4d023a8f 830 outVtx->SetIndices(nIndices,outindices);
44a76311 831 if (outindices) delete outindices;
fee9ef0d 832
919e537f 833 /*
834 printf("Vertex before removing tracks:");
fee9ef0d 835 inVtx->PrintStatus();
836 inVtx->PrintIndices();
919e537f 837 printf("Vertex after removing tracks:");
fee9ef0d 838 outVtx->PrintStatus();
839 outVtx->PrintIndices();
919e537f 840 */
fee9ef0d 841
842 return outVtx;
843}
844//---------------------------------------------------------------------------
6b908072 845AliESDVertex* AliVertexerTracks::RemoveConstraintFromVertex(AliESDVertex *inVtx,
846 Float_t *diamondxyz,
847 Float_t *diamondcov) const
848{
849//
850// Removes diamond constraint from fit of inVtx
851//
852
853 if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraint")) {
854 printf("ERROR: primary vertex has no constraint: cannot remove it\n");
855 return 0x0;
856 }
857 if(inVtx->GetNContributors()<3) {
858 printf("ERROR: primary vertex has less than 2 tracks: cannot remove contraint\n");
859 return 0x0;
860 }
861
862 // diamond constraint
863 TMatrixD vVb(3,3);
864 vVb(0,0) = diamondcov[0];
865 vVb(0,1) = diamondcov[1];
866 vVb(0,2) = 0.;
867 vVb(1,0) = diamondcov[1];
868 vVb(1,1) = diamondcov[2];
869 vVb(1,2) = 0.;
870 vVb(2,0) = 0.;
871 vVb(2,1) = 0.;
872 vVb(2,2) = diamondcov[5];
873 TMatrixD vVbInv(TMatrixD::kInverted,vVb);
874 TMatrixD rb(3,1);
875 rb(0,0) = diamondxyz[0];
876 rb(1,0) = diamondxyz[1];
877 rb(2,0) = diamondxyz[2];
878 TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb);
879
880 // input vertex
881 TMatrixD rv(3,1);
882 rv(0,0) = inVtx->GetXv();
883 rv(1,0) = inVtx->GetYv();
884 rv(2,0) = inVtx->GetZv();
885 TMatrixD vV(3,3);
886 Double_t cov[6];
887 inVtx->GetCovMatrix(cov);
888 vV(0,0) = cov[0];
889 vV(0,1) = cov[1]; vV(1,0) = cov[1];
890 vV(1,1) = cov[2];
891 vV(0,2) = cov[3]; vV(2,0) = cov[3];
892 vV(1,2) = cov[4]; vV(2,1) = cov[4];
893 vV(2,2) = cov[5];
894 TMatrixD vVInv(TMatrixD::kInverted,vV);
895 TMatrixD vVInvrv(vVInv,TMatrixD::kMult,rv);
896
897
898 TMatrixD sumWi = vVInv - vVbInv;
899
900
901 TMatrixD sumWiri = vVInvrv - vVbInvrb;
902
903 TMatrixD rvnew(3,1);
904 TMatrixD vVnew(3,3);
905
906 // new inverted of weights matrix
907 TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
908 vVnew = invsumWi;
909 // new position of primary vertex
910 rvnew.Mult(vVnew,sumWiri);
911
912 Double_t position[3];
913 position[0] = rvnew(0,0);
914 position[1] = rvnew(1,0);
915 position[2] = rvnew(2,0);
916 cov[0] = vVnew(0,0);
917 cov[1] = vVnew(0,1);
918 cov[2] = vVnew(1,1);
919 cov[3] = vVnew(0,2);
920 cov[4] = vVnew(1,2);
921 cov[5] = vVnew(2,2);
922
923
924 Double_t chi2 = inVtx->GetChi2();
925
926 // diamond constribution to chi2
927 TMatrixD deltar = rv; deltar -= rb;
928 TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar);
929 Double_t chi2b = deltar(0,0)*vVbInvdeltar(0,0)+
930 deltar(1,0)*vVbInvdeltar(1,0)+
931 deltar(2,0)*vVbInvdeltar(2,0);
932 // remove from total chi2
933 chi2 -= chi2b;
934
935 // store data in the vertex object
936 AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,inVtx->GetNContributors()-1);
937 outVtx->SetTitle("VertexerTracksNoConstraint");
938 UShort_t *inindices = inVtx->GetIndices();
939 Int_t nIndices = inVtx->GetNIndices();
940 outVtx->SetIndices(nIndices,inindices);
941
942 return outVtx;
943}
944//---------------------------------------------------------------------------
a00021a7 945void AliVertexerTracks::SetCuts(Double_t *cuts)
946{
947//
948// Cut values
949//
950 SetDCAcut(cuts[0]);
951 SetDCAcutIter0(cuts[1]);
952 SetMaxd0z0(cuts[2]);
5b78b791 953 if(fMode==0 && cuts[3]<0) SetITSrefitNotRequired();
954 SetMinClusters((Int_t)(TMath::Abs(cuts[3])));
a00021a7 955 SetMinTracks((Int_t)(cuts[4]));
956 SetNSigmad0(cuts[5]);
957 SetMinDetFitter(cuts[6]);
958 SetMaxTgl(cuts[7]);
959 SetFiducialRZ(cuts[8],cuts[9]);
8c75f668 960 fAlgo=(Int_t)(cuts[10]);
961 fAlgoIter0=(Int_t)(cuts[11]);
a00021a7 962
963 return;
964}
965//---------------------------------------------------------------------------
f09c879d 966void AliVertexerTracks::SetITSMode(Double_t dcacut,
967 Double_t dcacutIter0,
968 Double_t maxd0z0,
6b29d399 969 Int_t minCls,
f09c879d 970 Int_t mintrks,
971 Double_t nsigma,
972 Double_t mindetfitter,
6b29d399 973 Double_t maxtgl,
974 Double_t fidR,
8c75f668 975 Double_t fidZ,
976 Int_t finderAlgo,
977 Int_t finderAlgoIter0)
f09c879d 978{
979//
980// Cut values for ITS mode
981//
982 fMode = 0;
5b78b791 983 if(minCls>0) {
984 SetITSrefitRequired();
985 } else {
986 SetITSrefitNotRequired();
987 }
f09c879d 988 SetDCAcut(dcacut);
989 SetDCAcutIter0(dcacutIter0);
990 SetMaxd0z0(maxd0z0);
5b78b791 991 SetMinClusters(TMath::Abs(minCls));
f09c879d 992 SetMinTracks(mintrks);
993 SetNSigmad0(nsigma);
994 SetMinDetFitter(mindetfitter);
995 SetMaxTgl(maxtgl);
6b29d399 996 SetFiducialRZ(fidR,fidZ);
8c75f668 997 fAlgo=finderAlgo;
998 fAlgoIter0=finderAlgoIter0;
f09c879d 999
1000 return;
1001}
1002//---------------------------------------------------------------------------
1003void AliVertexerTracks::SetTPCMode(Double_t dcacut,
1004 Double_t dcacutIter0,
1005 Double_t maxd0z0,
6b29d399 1006 Int_t minCls,
f09c879d 1007 Int_t mintrks,
1008 Double_t nsigma,
1009 Double_t mindetfitter,
6b29d399 1010 Double_t maxtgl,
1011 Double_t fidR,
8c75f668 1012 Double_t fidZ,
1013 Int_t finderAlgo,
1014 Int_t finderAlgoIter0)
f09c879d 1015{
1016//
1017// Cut values for TPC mode
1018//
1019 fMode = 1;
1020 SetITSrefitNotRequired();
1021 SetDCAcut(dcacut);
1022 SetDCAcutIter0(dcacutIter0);
1023 SetMaxd0z0(maxd0z0);
6b29d399 1024 SetMinClusters(minCls);
f09c879d 1025 SetMinTracks(mintrks);
1026 SetNSigmad0(nsigma);
1027 SetMinDetFitter(mindetfitter);
1028 SetMaxTgl(maxtgl);
6b29d399 1029 SetFiducialRZ(fidR,fidZ);
8c75f668 1030 fAlgo=finderAlgo;
1031 fAlgoIter0=finderAlgoIter0;
f09c879d 1032
1033 return;
1034}
1035//---------------------------------------------------------------------------
fee9ef0d 1036void AliVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped)
1037{
1038//
1039// Mark the tracks not to be used in the vertex reconstruction.
1040// Tracks are identified by AliESDtrack::GetID()
1041//
1042 fNTrksToSkip = n; fTrksToSkip = new Int_t[n];
50ff0bcd 1043 for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
1044 return;
1045}
1046//---------------------------------------------------------------------------
fee9ef0d 1047void AliVertexerTracks::SetVtxStart(AliESDVertex *vtx)
1048{
50ff0bcd 1049//
1050// Set initial vertex knowledge
1051//
1052 vtx->GetXYZ(fNominalPos);
1053 vtx->GetCovMatrix(fNominalCov);
fee9ef0d 1054 SetConstraintOn();
50ff0bcd 1055 return;
1056}
2d57349e 1057//---------------------------------------------------------------------------
fee9ef0d 1058void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
dc3719f3 1059{
6d8df534 1060 AliExternalTrackParam *track1;
6d8df534 1061 const Int_t knacc = (Int_t)fTrkArraySel.GetEntriesFast();
8e104736 1062 AliStrLine **linarray = new AliStrLine* [knacc];
dc3719f3 1063 for(Int_t i=0; i<knacc; i++){
6d8df534 1064 track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
dc3719f3 1065 Double_t alpha=track1->GetAlpha();
1066 Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
1067 Double_t pos[3],dir[3],sigmasq[3];
f09c879d 1068 track1->GetXYZAt(mindist,GetFieldkG(),pos);
1069 track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
dc3719f3 1070 sigmasq[0]=TMath::Sin(alpha)*TMath::Sin(alpha)*track1->GetSigmaY2();
1071 sigmasq[1]=TMath::Cos(alpha)*TMath::Cos(alpha)*track1->GetSigmaY2();
1072 sigmasq[2]=track1->GetSigmaZ2();
146c29df 1073 TMatrixD ri(3,1);
1074 TMatrixD wWi(3,3);
8c75f668 1075 if(!TrackToPoint(track1,ri,wWi)) {optUseWeights=kFALSE;printf("WARNING\n");}
146c29df 1076 Double_t wmat[9];
1077 Int_t iel=0;
1078 for(Int_t ia=0;ia<3;ia++){
1079 for(Int_t ib=0;ib<3;ib++){
1080 wmat[iel]=wWi(ia,ib);
1081 iel++;
1082 }
1083 }
8e104736 1084 linarray[i] = new AliStrLine(pos,sigmasq,wmat,dir);
dc3719f3 1085 }
8e104736 1086 fVert=TrackletVertexFinder(linarray,knacc,optUseWeights);
1087 for(Int_t i=0; i<knacc; i++) delete linarray[i];
1088 delete [] linarray;
dc3719f3 1089}
1090//---------------------------------------------------------------------------
146c29df 1091AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t optUseWeights)
fee9ef0d 1092{
8e104736 1093 // Calculate the point at minimum distance to prepared tracks (TClonesArray)
dc3719f3 1094 const Int_t knacc = (Int_t)lines->GetEntriesFast();
8e104736 1095 AliStrLine** lines2 = new AliStrLine* [knacc];
1096 for(Int_t i=0; i<knacc; i++){
1097 lines2[i]= (AliStrLine*)lines->At(i);
1098 }
1099 AliESDVertex vert = TrackletVertexFinder(lines2,knacc,optUseWeights);
1100 delete [] lines2;
1101 return vert;
1102}
1103
1104//---------------------------------------------------------------------------
1105AliESDVertex AliVertexerTracks::TrackletVertexFinder(AliStrLine **lines, const Int_t knacc, Int_t optUseWeights)
1106{
1107 // Calculate the point at minimum distance to prepared tracks (array of AliStrLine)
1108
dc3719f3 1109 Double_t initPos[3]={0.,0.,0.};
2d57349e 1110
2d57349e 1111 Double_t (*vectP0)[3]=new Double_t [knacc][3];
1112 Double_t (*vectP1)[3]=new Double_t [knacc][3];
1113
1114 Double_t sum[3][3];
1115 Double_t dsum[3]={0,0,0};
146c29df 1116 TMatrixD sumWi(3,3);
1117 for(Int_t i=0;i<3;i++){
1118 for(Int_t j=0;j<3;j++){
1119 sum[i][j]=0;
1120 sumWi(i,j)=0.;
1121 }
1122 }
1123
2d57349e 1124 for(Int_t i=0; i<knacc; i++){
8e104736 1125 AliStrLine *line1 = lines[i];
dc3719f3 1126 Double_t p0[3],cd[3],sigmasq[3];
146c29df 1127 Double_t wmat[9];
8c75f668 1128 if(!line1) printf("ERROR %d %d\n",i,knacc);
2d57349e 1129 line1->GetP0(p0);
1130 line1->GetCd(cd);
dc3719f3 1131 line1->GetSigma2P0(sigmasq);
146c29df 1132 line1->GetWMatrix(wmat);
1133 TMatrixD wWi(3,3);
1134 Int_t iel=0;
1135 for(Int_t ia=0;ia<3;ia++){
1136 for(Int_t ib=0;ib<3;ib++){
1137 wWi(ia,ib)=wmat[iel];
1138 iel++;
1139 }
1140 }
1141
1142 sumWi+=wWi;
1143
2d57349e 1144 Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
1145 vectP0[i][0]=p0[0];
1146 vectP0[i][1]=p0[1];
1147 vectP0[i][2]=p0[2];
1148 vectP1[i][0]=p1[0];
1149 vectP1[i][1]=p1[1];
1150 vectP1[i][2]=p1[2];
1151
1152 Double_t matr[3][3];
1153 Double_t dknow[3];
1154 if(optUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
dc3719f3 1155 else GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
1156
2d57349e 1157
1158 for(Int_t iii=0;iii<3;iii++){
1159 dsum[iii]+=dknow[iii];
1160 for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
1161 }
2d57349e 1162 }
146c29df 1163
1164 TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1165 Double_t covmatrix[6];
1166 covmatrix[0] = invsumWi(0,0);
1167 covmatrix[1] = invsumWi(0,1);
1168 covmatrix[2] = invsumWi(1,1);
1169 covmatrix[3] = invsumWi(0,2);
1170 covmatrix[4] = invsumWi(1,2);
1171 covmatrix[5] = invsumWi(2,2);
1172
2d57349e 1173 Double_t vett[3][3];
1174 Double_t det=GetDeterminant3X3(sum);
dc3719f3 1175 Double_t sigma=0;
2d57349e 1176
fee9ef0d 1177 if(det!=0){
1178 for(Int_t zz=0;zz<3;zz++){
1179 for(Int_t ww=0;ww<3;ww++){
1180 for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
1181 }
1182 for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
1183 initPos[zz]=GetDeterminant3X3(vett)/det;
1184 }
1185
1186
1187 for(Int_t i=0; i<knacc; i++){
1188 Double_t p0[3]={0,0,0},p1[3]={0,0,0};
1189 for(Int_t ii=0;ii<3;ii++){
1190 p0[ii]=vectP0[i][ii];
1191 p1[ii]=vectP1[i][ii];
1192 }
1193 sigma+=GetStrLinMinDist(p0,p1,initPos);
1194 }
1195
f09c879d 1196 if(sigma>0.) {sigma=TMath::Sqrt(sigma);}else{sigma=999;}
fee9ef0d 1197 }else{
50ff0bcd 1198 sigma=999;
1199 }
146c29df 1200 AliESDVertex theVert(initPos,covmatrix,99999.,knacc);
1201 theVert.SetDispersion(sigma);
50ff0bcd 1202 delete vectP0;
1203 delete vectP1;
dc3719f3 1204 return theVert;
50ff0bcd 1205}
8e104736 1206
50ff0bcd 1207//---------------------------------------------------------------------------
6d8df534 1208Bool_t AliVertexerTracks::TrackToPoint(AliExternalTrackParam *t,
0bf2e2b4 1209 TMatrixD &ri,TMatrixD &wWi,
1210 Bool_t uUi3by3) const
fee9ef0d 1211{
1212//
6d8df534 1213// Extract from the AliExternalTrackParam the global coordinates ri and covariance matrix
fee9ef0d 1214// wWi of the space point that it represents (to be used in VertexFitter())
1215//
1216
1217
1218 Double_t rotAngle = t->GetAlpha();
1219 if(rotAngle<0.) rotAngle += 2.*TMath::Pi();
1220 Double_t cosRot = TMath::Cos(rotAngle);
1221 Double_t sinRot = TMath::Sin(rotAngle);
1222
1223 ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
1224 ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
1225 ri(2,0) = t->GetZ();
1226
0bf2e2b4 1227
1228
1229 if(!uUi3by3) {
1230 // matrix to go from global (x,y,z) to local (y,z);
1231 TMatrixD qQi(2,3);
1232 qQi(0,0) = -sinRot;
1233 qQi(0,1) = cosRot;
1234 qQi(0,2) = 0.;
1235 qQi(1,0) = 0.;
1236 qQi(1,1) = 0.;
1237 qQi(1,2) = 1.;
1238
1239 // covariance matrix of local (y,z) - inverted
0bf2e2b4 1240 TMatrixD uUi(2,2);
6d8df534 1241 uUi(0,0) = t->GetSigmaY2();
1242 uUi(0,1) = t->GetSigmaZY();
1243 uUi(1,0) = t->GetSigmaZY();
1244 uUi(1,1) = t->GetSigmaZ2();
919e537f 1245 //printf(" Ui :");
1246 //printf(" %f %f",uUi(0,0),uUi(0,1));
1247 //printf(" %f %f",uUi(1,0),uUi(1,1));
0bf2e2b4 1248
1249 if(uUi.Determinant() <= 0.) return kFALSE;
1250 TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1251
1252 // weights matrix: wWi = qQiT * uUiInv * qQi
1253 TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1254 TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1255 wWi = m;
1256 } else {
1257 if(fVert.GetNContributors()<1) AliFatal("Vertex from finder is empty");
1258 // matrix to go from global (x,y,z) to local (x,y,z);
1259 TMatrixD qQi(3,3);
1260 qQi(0,0) = cosRot;
1261 qQi(0,1) = sinRot;
1262 qQi(0,2) = 0.;
1263 qQi(1,0) = -sinRot;
1264 qQi(1,1) = cosRot;
1265 qQi(1,2) = 0.;
1266 qQi(2,0) = 0.;
1267 qQi(2,1) = 0.;
1268 qQi(2,2) = 1.;
1269
1270 // covariance of fVert along the track
1271 Double_t p[3],pt,ptot;
1272 t->GetPxPyPz(p);
1273 pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
1274 ptot = TMath::Sqrt(pt*pt+p[2]*p[2]);
3b768a1d 1275 Double_t cphi = p[0]/pt; //cos(phi)=px/pt
1276 Double_t sphi = p[1]/pt; //sin(phi)=py/pt
1277 Double_t clambda = pt/ptot; //cos(lambda)=pt/ptot
1278 Double_t slambda = p[2]/ptot; //sin(lambda)=pz/ptot
0bf2e2b4 1279 Double_t covfVert[6];
1280 fVert.GetCovMatrix(covfVert);
1281 Double_t covfVertalongt =
3b768a1d 1282 covfVert[0]*cphi*cphi*clambda*clambda
1283 +covfVert[1]*2.*cphi*sphi*clambda*clambda
1284 +covfVert[3]*2.*cphi*clambda*slambda
1285 +covfVert[2]*sphi*sphi*clambda*clambda
1286 +covfVert[4]*2.*sphi*clambda*slambda
1287 +covfVert[5]*slambda*slambda;
0bf2e2b4 1288 // covariance matrix of local (x,y,z) - inverted
1289 TMatrixD uUi(3,3);
3b768a1d 1290 uUi(0,0) = covfVertalongt * fnSigmaForUi00 * fnSigmaForUi00;
919e537f 1291 AliDebug(1,Form("=====> sqrtUi00 cm %f",TMath::Sqrt(uUi(0,0))));
0bf2e2b4 1292 uUi(0,1) = 0.;
1293 uUi(0,2) = 0.;
1294 uUi(1,0) = 0.;
6d8df534 1295 uUi(1,1) = t->GetSigmaY2();
1296 uUi(1,2) = t->GetSigmaZY();
0bf2e2b4 1297 uUi(2,0) = 0.;
6d8df534 1298 uUi(2,1) = t->GetSigmaZY();
1299 uUi(2,2) = t->GetSigmaZ2();
0bf2e2b4 1300 //printf(" Ui :\n");
1301 //printf(" %f %f\n",uUi(0,0),uUi(0,1));
1302 //printf(" %f %f\n",uUi(1,0),uUi(1,1));
fee9ef0d 1303
0bf2e2b4 1304 if(uUi.Determinant() <= 0.) return kFALSE;
1305 TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1306
1307 // weights matrix: wWi = qQiT * uUiInv * qQi
1308 TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1309 TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1310 wWi = m;
1311 }
1312
fee9ef0d 1313
1314 return kTRUE;
1315}
1316//---------------------------------------------------------------------------
6d8df534 1317void AliVertexerTracks::TooFewTracks()
fee9ef0d 1318{
50ff0bcd 1319//
6d8df534 1320// When the number of tracks is < fMinTracks,
1321// deal with vertices not found and prepare to exit
50ff0bcd 1322//
919e537f 1323 AliDebug(1,"TooFewTracks");
2d57349e 1324
50ff0bcd 1325 Double_t pos[3],err[3];
50ff0bcd 1326 pos[0] = fNominalPos[0];
1327 err[0] = TMath::Sqrt(fNominalCov[0]);
1328 pos[1] = fNominalPos[1];
1329 err[1] = TMath::Sqrt(fNominalCov[2]);
6d8df534 1330 pos[2] = fNominalPos[2];
1331 err[2] = TMath::Sqrt(fNominalCov[5]);
1332 Int_t ncontr = (err[0]>1. ? -1 : -3);
5b78b791 1333 if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
50ff0bcd 1334 fCurrentVertex = new AliESDVertex(pos,err);
1335 fCurrentVertex->SetNContributors(ncontr);
2d57349e 1336
fee9ef0d 1337 if(fConstraint) {
1338 fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
1339 } else {
50ff0bcd 1340 fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
fee9ef0d 1341 }
2d57349e 1342
6d8df534 1343 if(!fTrkArraySel.IsEmpty()) fTrkArraySel.Delete();
44a76311 1344 if(fIdSel) {delete fIdSel; fIdSel=NULL;}
1345 if(fTrksToSkip) {delete fTrksToSkip; fTrksToSkip=NULL;}
6d8df534 1346
50ff0bcd 1347 return;
1348}
1349//---------------------------------------------------------------------------
fee9ef0d 1350void AliVertexerTracks::VertexFinder(Int_t optUseWeights)
1351{
2d57349e 1352
50ff0bcd 1353 // Get estimate of vertex position in (x,y) from tracks DCA
1354
1355 Double_t initPos[3];
1356 initPos[2] = 0.;
1357 for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
6d8df534 1358 Int_t nacc = (Int_t)fTrkArraySel.GetEntriesFast();
50ff0bcd 1359 Double_t aver[3]={0.,0.,0.};
1360 Double_t aversq[3]={0.,0.,0.};
1361 Double_t sigmasq[3]={0.,0.,0.};
1362 Double_t sigma=0;
1363 Int_t ncombi = 0;
6d8df534 1364 AliExternalTrackParam *track1;
1365 AliExternalTrackParam *track2;
50ff0bcd 1366 Double_t pos[3],dir[3];
1367 Double_t alpha,mindist;
2d57349e 1368
50ff0bcd 1369 for(Int_t i=0; i<nacc; i++){
6d8df534 1370 track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
50ff0bcd 1371 alpha=track1->GetAlpha();
1372 mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
f09c879d 1373 track1->GetXYZAt(mindist,GetFieldkG(),pos);
1374 track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
50ff0bcd 1375 AliStrLine *line1 = new AliStrLine(pos,dir);
2d57349e 1376
50ff0bcd 1377 // AliStrLine *line1 = new AliStrLine();
f09c879d 1378 // track1->ApproximateHelixWithLine(mindist,GetFieldkG(),line1);
50ff0bcd 1379
1380 for(Int_t j=i+1; j<nacc; j++){
6d8df534 1381 track2 = (AliExternalTrackParam*)fTrkArraySel.At(j);
50ff0bcd 1382 alpha=track2->GetAlpha();
1383 mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
f09c879d 1384 track2->GetXYZAt(mindist,GetFieldkG(),pos);
1385 track2->GetPxPyPzAt(mindist,GetFieldkG(),dir);
50ff0bcd 1386 AliStrLine *line2 = new AliStrLine(pos,dir);
1387 // AliStrLine *line2 = new AliStrLine();
f09c879d 1388 // track2->ApproximateHelixWithLine(mindist,GetFieldkG(),line2);
50ff0bcd 1389 Double_t distCA=line2->GetDCA(line1);
fee9ef0d 1390 //printf("%d %d %f\n",i,j,distCA);
1391 if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
50ff0bcd 1392 Double_t pnt1[3],pnt2[3],crosspoint[3];
2d57349e 1393
50ff0bcd 1394 if(optUseWeights<=0){
1395 Int_t retcode = line2->Cross(line1,crosspoint);
1396 if(retcode>=0){
1397 ncombi++;
1398 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1399 for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1400 }
1401 }
1402 if(optUseWeights>0){
1403 Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2);
1404 if(retcode>=0){
5c03260f 1405 Double_t cs, sn;
50ff0bcd 1406 alpha=track1->GetAlpha();
1407 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1408 Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
1409 alpha=track2->GetAlpha();
1410 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1411 Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
1412 Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
1413 Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
1414 Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
1415 Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
1416 crosspoint[0]=wx1*pnt1[0] + wx2*pnt2[0];
1417 crosspoint[1]=wy1*pnt1[1] + wy2*pnt2[1];
1418 crosspoint[2]=wz1*pnt1[2] + wz2*pnt2[2];
1419
1420 ncombi++;
1421 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1422 for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1423 }
1424 }
1425 }
1426 delete line2;
1427 }
1428 delete line1;
71d84967 1429 }
50ff0bcd 1430 if(ncombi>0){
1431 for(Int_t jj=0;jj<3;jj++){
1432 initPos[jj] = aver[jj]/ncombi;
fee9ef0d 1433 //printf("%f\n",initPos[jj]);
50ff0bcd 1434 aversq[jj]/=ncombi;
1435 sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
1436 sigma+=sigmasq[jj];
1437 }
1438 sigma=TMath::Sqrt(TMath::Abs(sigma));
1439 }
1440 else {
1441 Warning("VertexFinder","Finder did not succed");
1442 sigma=999;
1443 }
1444 fVert.SetXYZ(initPos);
1445 fVert.SetDispersion(sigma);
1446 fVert.SetNContributors(ncombi);
2d57349e 1447}
ec1be5d5 1448//---------------------------------------------------------------------------
6d8df534 1449void AliVertexerTracks::VertexFitter()
fee9ef0d 1450{
ec1be5d5 1451//
1452// The optimal estimate of the vertex position is given by a "weighted
0bf2e2b4 1453// average of tracks positions".
fee9ef0d 1454// Original method: V. Karimaki, CMS Note 97/0051
ec1be5d5 1455//
6d8df534 1456 Bool_t useConstraint = fConstraint;
ec1be5d5 1457 Double_t initPos[3];
0bf2e2b4 1458 if(!fOnlyFitter) {
1459 fVert.GetXYZ(initPos);
1460 } else {
1461 initPos[0]=fNominalPos[0];
1462 initPos[1]=fNominalPos[1];
1463 initPos[2]=fNominalPos[2];
1464 }
1465
6d8df534 1466 Int_t nTrksSel = (Int_t)fTrkArraySel.GetEntries();
1467 if(nTrksSel==1) useConstraint=kTRUE;
919e537f 1468 AliDebug(1,"--- VertexFitter(): start");
1469 AliDebug(1,Form(" Number of tracks in array: %d\n",nTrksSel));
1470 AliDebug(1,Form(" Minimum # tracks required in fit: %d\n",fMinTracks));
1471 AliDebug(1,Form(" Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]));
1472 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 1473
0bf2e2b4 1474 // special treatment for few-tracks fits (e.g. secondary vertices)
6d8df534 1475 Bool_t uUi3by3 = kFALSE; if(nTrksSel<5 && !useConstraint) uUi3by3 = kTRUE;
0bf2e2b4 1476
ec1be5d5 1477 Int_t i,j,k,step=0;
1478 TMatrixD rv(3,1);
1479 TMatrixD vV(3,3);
1480 rv(0,0) = initPos[0];
1481 rv(1,0) = initPos[1];
1482 rv(2,0) = 0.;
1483 Double_t xlStart,alpha;
6d8df534 1484 Int_t nTrksUsed;
fee9ef0d 1485 Double_t chi2,chi2i,chi2b;
6d8df534 1486 AliExternalTrackParam *t = 0;
ec1be5d5 1487 Int_t failed = 0;
1488
50ff0bcd 1489 // initial vertex covariance matrix
1490 TMatrixD vVb(3,3);
1491 vVb(0,0) = fNominalCov[0];
1492 vVb(0,1) = fNominalCov[1];
1493 vVb(0,2) = 0.;
1494 vVb(1,0) = fNominalCov[1];
1495 vVb(1,1) = fNominalCov[2];
1496 vVb(1,2) = 0.;
1497 vVb(2,0) = 0.;
1498 vVb(2,1) = 0.;
1499 vVb(2,2) = fNominalCov[5];
1500 TMatrixD vVbInv(TMatrixD::kInverted,vVb);
1501 TMatrixD rb(3,1);
1502 rb(0,0) = fNominalPos[0];
1503 rb(1,0) = fNominalPos[1];
1504 rb(2,0) = fNominalPos[2];
1505 TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb);
1506
1507
eae1677e 1508 // 2 steps:
1509 // 1st - estimate of vtx using all tracks
1510 // 2nd - estimate of global chi2
1511 for(step=0; step<2; step++) {
919e537f 1512 AliDebug(1,Form(" step = %d\n",step));
ec1be5d5 1513 chi2 = 0.;
6d8df534 1514 nTrksUsed = 0;
ec1be5d5 1515
0bf2e2b4 1516 if(step==1) { initPos[0]=rv(0,0); initPos[0]=rv(1,0); }
1517
ec1be5d5 1518 TMatrixD sumWiri(3,1);
1519 TMatrixD sumWi(3,3);
1520 for(i=0; i<3; i++) {
1521 sumWiri(i,0) = 0.;
1522 for(j=0; j<3; j++) sumWi(j,i) = 0.;
1523 }
1524
fee9ef0d 1525 // mean vertex constraint
1526 if(useConstraint) {
50ff0bcd 1527 for(i=0;i<3;i++) {
1528 sumWiri(i,0) += vVbInvrb(i,0);
1529 for(k=0;k<3;k++) sumWi(i,k) += vVbInv(i,k);
07680cae 1530 }
fee9ef0d 1531 // chi2
1532 TMatrixD deltar = rv; deltar -= rb;
1533 TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar);
1534 chi2b = deltar(0,0)*vVbInvdeltar(0,0)+
1535 deltar(1,0)*vVbInvdeltar(1,0)+
1536 deltar(2,0)*vVbInvdeltar(2,0);
1537 chi2 += chi2b;
07680cae 1538 }
1539
ec1be5d5 1540 // loop on tracks
6d8df534 1541 for(k=0; k<nTrksSel; k++) {
ec1be5d5 1542 // get track from track array
6d8df534 1543 t = (AliExternalTrackParam*)fTrkArraySel.At(k);
ec1be5d5 1544 alpha = t->GetAlpha();
1545 xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
fee9ef0d 1546 // to vtxSeed (from finder)
6d8df534 1547 t->PropagateTo(xlStart,GetFieldkG());
fee9ef0d 1548
ec1be5d5 1549 // vector of track global coordinates
1550 TMatrixD ri(3,1);
fee9ef0d 1551 // covariance matrix of ri
1552 TMatrixD wWi(3,3);
1553
1554 // get space point from track
0bf2e2b4 1555 if(!TrackToPoint(t,ri,wWi,uUi3by3)) continue;
fee9ef0d 1556 TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
ec1be5d5 1557
1558 // track chi2
1559 TMatrixD deltar = rv; deltar -= ri;
1560 TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
1561 chi2i = deltar(0,0)*wWideltar(0,0)+
1562 deltar(1,0)*wWideltar(1,0)+
1563 deltar(2,0)*wWideltar(2,0);
1564
ec1be5d5 1565 // add to total chi2
1566 chi2 += chi2i;
1567
ec1be5d5 1568 sumWiri += wWiri;
1569 sumWi += wWi;
1570
6d8df534 1571 nTrksUsed++;
ec1be5d5 1572 } // end loop on tracks
1573
6d8df534 1574 if(nTrksUsed < fMinTracks) {
ec1be5d5 1575 failed=1;
1576 continue;
1577 }
1578
1579 Double_t determinant = sumWi.Determinant();
f09c879d 1580 if(determinant < fMinDetFitter) {
919e537f 1581 AliDebug(1,Form("det(V) = %f (<%f)\n",determinant,fMinDetFitter));
ec1be5d5 1582 failed=1;
1583 continue;
1584 }
1585
0bf2e2b4 1586 if(step==0) {
1587 // inverted of weights matrix
1588 TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1589 vV = invsumWi;
1590 // position of primary vertex
1591 rv.Mult(vV,sumWiri);
1592 }
eae1677e 1593 } // end loop on the 2 steps
ec1be5d5 1594
ec1be5d5 1595 if(failed) {
6d8df534 1596 TooFewTracks();
ec1be5d5 1597 return;
1598 }
1599
1600 Double_t position[3];
1601 position[0] = rv(0,0);
1602 position[1] = rv(1,0);
1603 position[2] = rv(2,0);
1604 Double_t covmatrix[6];
1605 covmatrix[0] = vV(0,0);
1606 covmatrix[1] = vV(0,1);
1607 covmatrix[2] = vV(1,1);
1608 covmatrix[3] = vV(0,2);
1609 covmatrix[4] = vV(1,2);
1610 covmatrix[5] = vV(2,2);
1611
dc3719f3 1612 // for correct chi2/ndf, count constraint as additional "track"
6d8df534 1613 if(fConstraint) nTrksUsed++;
ec1be5d5 1614 // store data in the vertex object
5b78b791 1615 if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
6d8df534 1616 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
ec1be5d5 1617
919e537f 1618 AliDebug(1," Vertex after fit:");
1619 AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
1620 AliDebug(1,"--- VertexFitter(): finish\n");
1621
ec1be5d5 1622
1623 return;
1624}
50ff0bcd 1625//----------------------------------------------------------------------------
6d8df534 1626AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray,
1627 UShort_t *id,
1628 Bool_t optUseFitter,
97669588 1629 Bool_t optPropagate,
1630 Bool_t optUseDiamondConstraint)
fee9ef0d 1631{
eae1677e 1632//
6d8df534 1633// Return vertex from tracks (AliExternalTrackParam) in array
eae1677e 1634//
5b78b791 1635 fCurrentVertex = 0;
97669588 1636 // set optUseDiamondConstraint=TRUE only if you are reconstructing the
1637 // primary vertex!
1638 if(optUseDiamondConstraint) {
1639 SetConstraintOn();
1640 } else {
1641 SetConstraintOff();
1642 }
fee9ef0d 1643
50ff0bcd 1644 // get tracks and propagate them to initial vertex position
6d8df534 1645 fIdSel = new UShort_t[(Int_t)trkArray->GetEntriesFast()];
1646 Int_t nTrksSel = PrepareTracks(*trkArray,id,0);
97669588 1647 if((!optUseDiamondConstraint && nTrksSel<TMath::Max(2,fMinTracks)) ||
1648 (optUseDiamondConstraint && nTrksSel<1)) {
6d8df534 1649 TooFewTracks();
fee9ef0d 1650 return fCurrentVertex;
50ff0bcd 1651 }
1652
6d8df534 1653 // vertex finder
97669588 1654 if(nTrksSel==1) {
1655 AliDebug(1,"Just one track");
1656 OneTrackVertFinder();
1657 } else {
1658 switch (fAlgo) {
1659 case 1: StrLinVertexFinderMinDist(1); break;
1660 case 2: StrLinVertexFinderMinDist(0); break;
1661 case 3: HelixVertexFinder(); break;
1662 case 4: VertexFinder(1); break;
1663 case 5: VertexFinder(0); break;
1664 default: printf("Wrong algorithm\n"); break;
1665 }
fee9ef0d 1666 }
919e537f 1667 AliDebug(1," Vertex finding completed\n");
fee9ef0d 1668
1669 // vertex fitter
6d8df534 1670 if(optUseFitter) {
1671 VertexFitter();
1672 } else {
fee9ef0d 1673 Double_t position[3]={fVert.GetXv(),fVert.GetYv(),fVert.GetZv()};
146c29df 1674 Double_t covmatrix[6];
1675 fVert.GetCovMatrix(covmatrix);
fee9ef0d 1676 Double_t chi2=99999.;
6d8df534 1677 Int_t nTrksUsed=fVert.GetNContributors();
1678 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
fee9ef0d 1679 }
1680 fCurrentVertex->SetDispersion(fVert.GetDispersion());
1681
1682
1683 // set indices of used tracks and propagate track to found vertex
50ff0bcd 1684 UShort_t *indices = 0;
6d8df534 1685 Double_t d0z0[2],covd0z0[3];
1686 AliExternalTrackParam *t = 0;
fee9ef0d 1687 if(fCurrentVertex->GetNContributors()>0) {
44a76311 1688 indices = new UShort_t[fTrkArraySel.GetEntriesFast()];
6d8df534 1689 for(Int_t jj=0; jj<(Int_t)fTrkArraySel.GetEntriesFast(); jj++) {
1690 indices[jj] = fIdSel[jj];
1691 t = (AliExternalTrackParam*)fTrkArraySel.At(jj);
1692 if(optPropagate && optUseFitter) {
fee9ef0d 1693 if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
6d8df534 1694 t->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
919e537f 1695 AliDebug(1,Form("Track %d propagated to found vertex",jj));
6d8df534 1696 } else {
fee9ef0d 1697 AliWarning("Found vertex outside beam pipe!");
1698 }
1699 }
50ff0bcd 1700 }
fee9ef0d 1701 fCurrentVertex->SetIndices(fCurrentVertex->GetNContributors(),indices);
50ff0bcd 1702 }
6d8df534 1703
1704 // clean up
44a76311 1705 if (indices) {delete indices; indices=NULL;}
1706 delete fIdSel; fIdSel=NULL;
6d8df534 1707 fTrkArraySel.Delete();
50ff0bcd 1708
fee9ef0d 1709 return fCurrentVertex;
eae1677e 1710}
50ff0bcd 1711//----------------------------------------------------------------------------
6d8df534 1712AliESDVertex* AliVertexerTracks::VertexForSelectedESDTracks(TObjArray *trkArray,
97669588 1713 Bool_t optUseFitter,
1714 Bool_t optPropagate,
1715 Bool_t optUseDiamondConstraint)
1716
fee9ef0d 1717{
50ff0bcd 1718//
6d8df534 1719// Return vertex from array of ESD tracks
50ff0bcd 1720//
5834e9a1 1721
6d8df534 1722 Int_t nTrks = (Int_t)trkArray->GetEntriesFast();
1723 UShort_t *id = new UShort_t[nTrks];
1724
1725 AliESDtrack *esdt = 0;
50ff0bcd 1726 for(Int_t i=0; i<nTrks; i++){
6d8df534 1727 esdt = (AliESDtrack*)trkArray->At(i);
1728 id[i] = (UShort_t)esdt->GetID();
50ff0bcd 1729 }
1730
97669588 1731 VertexForSelectedTracks(trkArray,id,optUseFitter,optPropagate,optUseDiamondConstraint);
6d8df534 1732
44a76311 1733 delete id; id=NULL;
6d8df534 1734
1735 return fCurrentVertex;
50ff0bcd 1736}
1737//--------------------------------------------------------------------------