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