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