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