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