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