]> 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
correcting raw <-> local transformation according to the 0.5 pad offset applied in...
[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//---------------------------------------------------------------------------
a00021a7 992void AliVertexerTracks::SetCuts(Double_t *cuts)
993{
994//
995// Cut values
996//
997 SetDCAcut(cuts[0]);
998 SetDCAcutIter0(cuts[1]);
999 SetMaxd0z0(cuts[2]);
5b78b791 1000 if(fMode==0 && cuts[3]<0) SetITSrefitNotRequired();
1001 SetMinClusters((Int_t)(TMath::Abs(cuts[3])));
a00021a7 1002 SetMinTracks((Int_t)(cuts[4]));
1003 SetNSigmad0(cuts[5]);
1004 SetMinDetFitter(cuts[6]);
1005 SetMaxTgl(cuts[7]);
1006 SetFiducialRZ(cuts[8],cuts[9]);
8c75f668 1007 fAlgo=(Int_t)(cuts[10]);
1008 fAlgoIter0=(Int_t)(cuts[11]);
3f2db92f 1009 //
1010 if (cuts[12]>1.) SetMVTukey2(cuts[12]);
1011 if (cuts[13]>1.) SetMVSig2Ini(cuts[13]);
1012 if (cuts[14]>0.1) SetMVMaxSigma2(cuts[14]);
1013 if (cuts[15]>1e-5) SetMVMinSig2Red(cuts[15]);
1014 if (cuts[16]>1e-5) SetMVMinDst(cuts[16]);
1015 if (cuts[17]>0.5) SetMVScanStep(cuts[17]);
1016 SetMVMaxWghNtr(cuts[18]);
1017 SetMVFinalWBinary(cuts[19]>0);
1018 if (cuts[20]>20.) SetBCSpacing(int(cuts[20]));
1019 //
1020 if (fAlgo==kMultiVertexer) SetSelectOnTOFBunchCrossing(kTRUE,kTRUE);
1021 else SetSelectOnTOFBunchCrossing(kFALSE,kTRUE);
a00021a7 1022 return;
1023}
1024//---------------------------------------------------------------------------
f09c879d 1025void AliVertexerTracks::SetITSMode(Double_t dcacut,
1026 Double_t dcacutIter0,
1027 Double_t maxd0z0,
6b29d399 1028 Int_t minCls,
f09c879d 1029 Int_t mintrks,
1030 Double_t nsigma,
1031 Double_t mindetfitter,
6b29d399 1032 Double_t maxtgl,
1033 Double_t fidR,
8c75f668 1034 Double_t fidZ,
1035 Int_t finderAlgo,
1036 Int_t finderAlgoIter0)
f09c879d 1037{
1038//
1039// Cut values for ITS mode
1040//
1041 fMode = 0;
5b78b791 1042 if(minCls>0) {
1043 SetITSrefitRequired();
1044 } else {
1045 SetITSrefitNotRequired();
1046 }
f09c879d 1047 SetDCAcut(dcacut);
1048 SetDCAcutIter0(dcacutIter0);
1049 SetMaxd0z0(maxd0z0);
5b78b791 1050 SetMinClusters(TMath::Abs(minCls));
f09c879d 1051 SetMinTracks(mintrks);
1052 SetNSigmad0(nsigma);
1053 SetMinDetFitter(mindetfitter);
1054 SetMaxTgl(maxtgl);
6b29d399 1055 SetFiducialRZ(fidR,fidZ);
8c75f668 1056 fAlgo=finderAlgo;
1057 fAlgoIter0=finderAlgoIter0;
f09c879d 1058
1059 return;
1060}
1061//---------------------------------------------------------------------------
1062void AliVertexerTracks::SetTPCMode(Double_t dcacut,
1063 Double_t dcacutIter0,
1064 Double_t maxd0z0,
6b29d399 1065 Int_t minCls,
f09c879d 1066 Int_t mintrks,
1067 Double_t nsigma,
1068 Double_t mindetfitter,
6b29d399 1069 Double_t maxtgl,
1070 Double_t fidR,
8c75f668 1071 Double_t fidZ,
1072 Int_t finderAlgo,
1073 Int_t finderAlgoIter0)
f09c879d 1074{
1075//
1076// Cut values for TPC mode
1077//
1078 fMode = 1;
1079 SetITSrefitNotRequired();
1080 SetDCAcut(dcacut);
1081 SetDCAcutIter0(dcacutIter0);
1082 SetMaxd0z0(maxd0z0);
6b29d399 1083 SetMinClusters(minCls);
f09c879d 1084 SetMinTracks(mintrks);
1085 SetNSigmad0(nsigma);
1086 SetMinDetFitter(mindetfitter);
1087 SetMaxTgl(maxtgl);
6b29d399 1088 SetFiducialRZ(fidR,fidZ);
8c75f668 1089 fAlgo=finderAlgo;
1090 fAlgoIter0=finderAlgoIter0;
f09c879d 1091
1092 return;
1093}
1094//---------------------------------------------------------------------------
f5740e9a 1095void AliVertexerTracks::SetSkipTracks(Int_t n,const Int_t *skipped)
fee9ef0d 1096{
1097//
1098// Mark the tracks not to be used in the vertex reconstruction.
1099// Tracks are identified by AliESDtrack::GetID()
1100//
1101 fNTrksToSkip = n; fTrksToSkip = new Int_t[n];
50ff0bcd 1102 for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
1103 return;
1104}
1105//---------------------------------------------------------------------------
fee9ef0d 1106void AliVertexerTracks::SetVtxStart(AliESDVertex *vtx)
1107{
50ff0bcd 1108//
1109// Set initial vertex knowledge
1110//
1111 vtx->GetXYZ(fNominalPos);
1112 vtx->GetCovMatrix(fNominalCov);
fee9ef0d 1113 SetConstraintOn();
50ff0bcd 1114 return;
1115}
2d57349e 1116//---------------------------------------------------------------------------
fee9ef0d 1117void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
dc3719f3 1118{
6d8df534 1119 AliExternalTrackParam *track1;
6d8df534 1120 const Int_t knacc = (Int_t)fTrkArraySel.GetEntriesFast();
8e104736 1121 AliStrLine **linarray = new AliStrLine* [knacc];
dc3719f3 1122 for(Int_t i=0; i<knacc; i++){
6d8df534 1123 track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
dc3719f3 1124 Double_t alpha=track1->GetAlpha();
1125 Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
1126 Double_t pos[3],dir[3],sigmasq[3];
f09c879d 1127 track1->GetXYZAt(mindist,GetFieldkG(),pos);
1128 track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
dc3719f3 1129 sigmasq[0]=TMath::Sin(alpha)*TMath::Sin(alpha)*track1->GetSigmaY2();
1130 sigmasq[1]=TMath::Cos(alpha)*TMath::Cos(alpha)*track1->GetSigmaY2();
1131 sigmasq[2]=track1->GetSigmaZ2();
146c29df 1132 TMatrixD ri(3,1);
1133 TMatrixD wWi(3,3);
8c75f668 1134 if(!TrackToPoint(track1,ri,wWi)) {optUseWeights=kFALSE;printf("WARNING\n");}
146c29df 1135 Double_t wmat[9];
1136 Int_t iel=0;
1137 for(Int_t ia=0;ia<3;ia++){
1138 for(Int_t ib=0;ib<3;ib++){
1139 wmat[iel]=wWi(ia,ib);
1140 iel++;
1141 }
1142 }
8e104736 1143 linarray[i] = new AliStrLine(pos,sigmasq,wmat,dir);
dc3719f3 1144 }
8e104736 1145 fVert=TrackletVertexFinder(linarray,knacc,optUseWeights);
1146 for(Int_t i=0; i<knacc; i++) delete linarray[i];
1147 delete [] linarray;
dc3719f3 1148}
1149//---------------------------------------------------------------------------
f5740e9a 1150AliESDVertex AliVertexerTracks::TrackletVertexFinder(const TClonesArray *lines, Int_t optUseWeights)
fee9ef0d 1151{
8e104736 1152 // Calculate the point at minimum distance to prepared tracks (TClonesArray)
dc3719f3 1153 const Int_t knacc = (Int_t)lines->GetEntriesFast();
8e104736 1154 AliStrLine** lines2 = new AliStrLine* [knacc];
1155 for(Int_t i=0; i<knacc; i++){
1156 lines2[i]= (AliStrLine*)lines->At(i);
1157 }
1158 AliESDVertex vert = TrackletVertexFinder(lines2,knacc,optUseWeights);
1159 delete [] lines2;
1160 return vert;
1161}
1162
1163//---------------------------------------------------------------------------
1164AliESDVertex AliVertexerTracks::TrackletVertexFinder(AliStrLine **lines, const Int_t knacc, Int_t optUseWeights)
1165{
1166 // Calculate the point at minimum distance to prepared tracks (array of AliStrLine)
1167
dc3719f3 1168 Double_t initPos[3]={0.,0.,0.};
2d57349e 1169
3f2db92f 1170 Double_t (*vectP0)[3]=new Double_t [knacc][3];
1171 Double_t (*vectP1)[3]=new Double_t [knacc][3];
2d57349e 1172
1173 Double_t sum[3][3];
1174 Double_t dsum[3]={0,0,0};
146c29df 1175 TMatrixD sumWi(3,3);
1176 for(Int_t i=0;i<3;i++){
1177 for(Int_t j=0;j<3;j++){
1178 sum[i][j]=0;
1179 sumWi(i,j)=0.;
1180 }
1181 }
1182
2d57349e 1183 for(Int_t i=0; i<knacc; i++){
8e104736 1184 AliStrLine *line1 = lines[i];
dc3719f3 1185 Double_t p0[3],cd[3],sigmasq[3];
146c29df 1186 Double_t wmat[9];
3f2db92f 1187 if(!line1) printf("ERROR %d %d\n",i,knacc);
2d57349e 1188 line1->GetP0(p0);
1189 line1->GetCd(cd);
dc3719f3 1190 line1->GetSigma2P0(sigmasq);
146c29df 1191 line1->GetWMatrix(wmat);
1192 TMatrixD wWi(3,3);
1193 Int_t iel=0;
1194 for(Int_t ia=0;ia<3;ia++){
1195 for(Int_t ib=0;ib<3;ib++){
1196 wWi(ia,ib)=wmat[iel];
1197 iel++;
1198 }
1199 }
1200
1201 sumWi+=wWi;
1202
2d57349e 1203 Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
1204 vectP0[i][0]=p0[0];
1205 vectP0[i][1]=p0[1];
1206 vectP0[i][2]=p0[2];
1207 vectP1[i][0]=p1[0];
1208 vectP1[i][1]=p1[1];
1209 vectP1[i][2]=p1[2];
1210
1211 Double_t matr[3][3];
1212 Double_t dknow[3];
1213 if(optUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
dc3719f3 1214 else GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
1215
2d57349e 1216
1217 for(Int_t iii=0;iii<3;iii++){
1218 dsum[iii]+=dknow[iii];
1219 for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
1220 }
2d57349e 1221 }
146c29df 1222
1223 TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1224 Double_t covmatrix[6];
1225 covmatrix[0] = invsumWi(0,0);
1226 covmatrix[1] = invsumWi(0,1);
1227 covmatrix[2] = invsumWi(1,1);
1228 covmatrix[3] = invsumWi(0,2);
1229 covmatrix[4] = invsumWi(1,2);
1230 covmatrix[5] = invsumWi(2,2);
1231
2d57349e 1232 Double_t vett[3][3];
1233 Double_t det=GetDeterminant3X3(sum);
dc3719f3 1234 Double_t sigma=0;
2d57349e 1235
f5740e9a 1236 if(TMath::Abs(det) > kAlmost0){
fee9ef0d 1237 for(Int_t zz=0;zz<3;zz++){
1238 for(Int_t ww=0;ww<3;ww++){
1239 for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
1240 }
1241 for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
1242 initPos[zz]=GetDeterminant3X3(vett)/det;
1243 }
1244
1245
1246 for(Int_t i=0; i<knacc; i++){
1247 Double_t p0[3]={0,0,0},p1[3]={0,0,0};
1248 for(Int_t ii=0;ii<3;ii++){
1249 p0[ii]=vectP0[i][ii];
1250 p1[ii]=vectP1[i][ii];
1251 }
1252 sigma+=GetStrLinMinDist(p0,p1,initPos);
1253 }
1254
f09c879d 1255 if(sigma>0.) {sigma=TMath::Sqrt(sigma);}else{sigma=999;}
fee9ef0d 1256 }else{
50ff0bcd 1257 sigma=999;
1258 }
146c29df 1259 AliESDVertex theVert(initPos,covmatrix,99999.,knacc);
1260 theVert.SetDispersion(sigma);
4ce766eb 1261 delete [] vectP0;
1262 delete [] vectP1;
dc3719f3 1263 return theVert;
50ff0bcd 1264}
8e104736 1265
50ff0bcd 1266//---------------------------------------------------------------------------
6d8df534 1267Bool_t AliVertexerTracks::TrackToPoint(AliExternalTrackParam *t,
817e1004 1268 TMatrixD &ri,TMatrixD &wWi,
0bf2e2b4 1269 Bool_t uUi3by3) const
fee9ef0d 1270{
1271//
6d8df534 1272// Extract from the AliExternalTrackParam the global coordinates ri and covariance matrix
fee9ef0d 1273// wWi of the space point that it represents (to be used in VertexFitter())
1274//
1275
1276
1277 Double_t rotAngle = t->GetAlpha();
1278 if(rotAngle<0.) rotAngle += 2.*TMath::Pi();
1279 Double_t cosRot = TMath::Cos(rotAngle);
1280 Double_t sinRot = TMath::Sin(rotAngle);
3f2db92f 1281 /*
1282 // RS >>>
1283 Double_t lambda = TMath::ATan(t->GetTgl());
1284 Double_t cosLam = TMath::Cos(lambda);
1285 Double_t sinLam = TMath::Sin(lambda);
1286 // RS <<<
1287 */
fee9ef0d 1288 ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
1289 ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
1290 ri(2,0) = t->GetZ();
1291
0bf2e2b4 1292 if(!uUi3by3) {
1293 // matrix to go from global (x,y,z) to local (y,z);
1294 TMatrixD qQi(2,3);
1295 qQi(0,0) = -sinRot;
1296 qQi(0,1) = cosRot;
1297 qQi(0,2) = 0.;
3f2db92f 1298 //
0bf2e2b4 1299 qQi(1,0) = 0.;
1300 qQi(1,1) = 0.;
1301 qQi(1,2) = 1.;
3f2db92f 1302 //
1303 // RS: Added polar inclination
1304 /*
1305 qQi(1,0) = -sinLam*cosRot;
1306 qQi(1,1) = -sinLam*sinRot;
1307 qQi(1,2) = cosLam;
1308 */
0bf2e2b4 1309 // covariance matrix of local (y,z) - inverted
0bf2e2b4 1310 TMatrixD uUi(2,2);
6d8df534 1311 uUi(0,0) = t->GetSigmaY2();
1312 uUi(0,1) = t->GetSigmaZY();
1313 uUi(1,0) = t->GetSigmaZY();
1314 uUi(1,1) = t->GetSigmaZ2();
919e537f 1315 //printf(" Ui :");
1316 //printf(" %f %f",uUi(0,0),uUi(0,1));
1317 //printf(" %f %f",uUi(1,0),uUi(1,1));
0bf2e2b4 1318
1319 if(uUi.Determinant() <= 0.) return kFALSE;
1320 TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1321
1322 // weights matrix: wWi = qQiT * uUiInv * qQi
1323 TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1324 TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1325 wWi = m;
1326 } else {
1327 if(fVert.GetNContributors()<1) AliFatal("Vertex from finder is empty");
1328 // matrix to go from global (x,y,z) to local (x,y,z);
1329 TMatrixD qQi(3,3);
1330 qQi(0,0) = cosRot;
1331 qQi(0,1) = sinRot;
1332 qQi(0,2) = 0.;
1333 qQi(1,0) = -sinRot;
1334 qQi(1,1) = cosRot;
1335 qQi(1,2) = 0.;
1336 qQi(2,0) = 0.;
1337 qQi(2,1) = 0.;
1338 qQi(2,2) = 1.;
1339
1340 // covariance of fVert along the track
1341 Double_t p[3],pt,ptot;
1342 t->GetPxPyPz(p);
1343 pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
1344 ptot = TMath::Sqrt(pt*pt+p[2]*p[2]);
3b768a1d 1345 Double_t cphi = p[0]/pt; //cos(phi)=px/pt
1346 Double_t sphi = p[1]/pt; //sin(phi)=py/pt
1347 Double_t clambda = pt/ptot; //cos(lambda)=pt/ptot
1348 Double_t slambda = p[2]/ptot; //sin(lambda)=pz/ptot
0bf2e2b4 1349 Double_t covfVert[6];
1350 fVert.GetCovMatrix(covfVert);
1351 Double_t covfVertalongt =
3b768a1d 1352 covfVert[0]*cphi*cphi*clambda*clambda
1353 +covfVert[1]*2.*cphi*sphi*clambda*clambda
1354 +covfVert[3]*2.*cphi*clambda*slambda
1355 +covfVert[2]*sphi*sphi*clambda*clambda
1356 +covfVert[4]*2.*sphi*clambda*slambda
1357 +covfVert[5]*slambda*slambda;
0bf2e2b4 1358 // covariance matrix of local (x,y,z) - inverted
1359 TMatrixD uUi(3,3);
3b768a1d 1360 uUi(0,0) = covfVertalongt * fnSigmaForUi00 * fnSigmaForUi00;
919e537f 1361 AliDebug(1,Form("=====> sqrtUi00 cm %f",TMath::Sqrt(uUi(0,0))));
0bf2e2b4 1362 uUi(0,1) = 0.;
1363 uUi(0,2) = 0.;
1364 uUi(1,0) = 0.;
6d8df534 1365 uUi(1,1) = t->GetSigmaY2();
1366 uUi(1,2) = t->GetSigmaZY();
0bf2e2b4 1367 uUi(2,0) = 0.;
6d8df534 1368 uUi(2,1) = t->GetSigmaZY();
1369 uUi(2,2) = t->GetSigmaZ2();
0bf2e2b4 1370 //printf(" Ui :\n");
1371 //printf(" %f %f\n",uUi(0,0),uUi(0,1));
1372 //printf(" %f %f\n",uUi(1,0),uUi(1,1));
fee9ef0d 1373
0bf2e2b4 1374 if(uUi.Determinant() <= 0.) return kFALSE;
1375 TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1376
1377 // weights matrix: wWi = qQiT * uUiInv * qQi
1378 TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1379 TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1380 wWi = m;
1381 }
1382
fee9ef0d 1383
1384 return kTRUE;
1385}
1386//---------------------------------------------------------------------------
6d8df534 1387void AliVertexerTracks::TooFewTracks()
fee9ef0d 1388{
50ff0bcd 1389//
6d8df534 1390// When the number of tracks is < fMinTracks,
1391// deal with vertices not found and prepare to exit
50ff0bcd 1392//
919e537f 1393 AliDebug(1,"TooFewTracks");
2d57349e 1394
50ff0bcd 1395 Double_t pos[3],err[3];
50ff0bcd 1396 pos[0] = fNominalPos[0];
1397 err[0] = TMath::Sqrt(fNominalCov[0]);
1398 pos[1] = fNominalPos[1];
1399 err[1] = TMath::Sqrt(fNominalCov[2]);
6d8df534 1400 pos[2] = fNominalPos[2];
1401 err[2] = TMath::Sqrt(fNominalCov[5]);
1402 Int_t ncontr = (err[0]>1. ? -1 : -3);
5b78b791 1403 if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
50ff0bcd 1404 fCurrentVertex = new AliESDVertex(pos,err);
1405 fCurrentVertex->SetNContributors(ncontr);
2d57349e 1406
fee9ef0d 1407 if(fConstraint) {
1408 fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
1409 } else {
50ff0bcd 1410 fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
fee9ef0d 1411 }
2d57349e 1412
6d8df534 1413 if(!fTrkArraySel.IsEmpty()) fTrkArraySel.Delete();
6fefa5ce 1414 if(fIdSel) {delete [] fIdSel; fIdSel=NULL;}
1415 if(fTrksToSkip) {delete [] fTrksToSkip; fTrksToSkip=NULL;}
6d8df534 1416
50ff0bcd 1417 return;
1418}
1419//---------------------------------------------------------------------------
fee9ef0d 1420void AliVertexerTracks::VertexFinder(Int_t optUseWeights)
1421{
2d57349e 1422
50ff0bcd 1423 // Get estimate of vertex position in (x,y) from tracks DCA
1424
1425 Double_t initPos[3];
1426 initPos[2] = 0.;
1427 for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
6d8df534 1428 Int_t nacc = (Int_t)fTrkArraySel.GetEntriesFast();
50ff0bcd 1429 Double_t aver[3]={0.,0.,0.};
1430 Double_t aversq[3]={0.,0.,0.};
1431 Double_t sigmasq[3]={0.,0.,0.};
1432 Double_t sigma=0;
1433 Int_t ncombi = 0;
6d8df534 1434 AliExternalTrackParam *track1;
1435 AliExternalTrackParam *track2;
50ff0bcd 1436 Double_t pos[3],dir[3];
1437 Double_t alpha,mindist;
2d57349e 1438
50ff0bcd 1439 for(Int_t i=0; i<nacc; i++){
6d8df534 1440 track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
50ff0bcd 1441 alpha=track1->GetAlpha();
1442 mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
f09c879d 1443 track1->GetXYZAt(mindist,GetFieldkG(),pos);
1444 track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
50ff0bcd 1445 AliStrLine *line1 = new AliStrLine(pos,dir);
2d57349e 1446
50ff0bcd 1447 // AliStrLine *line1 = new AliStrLine();
f09c879d 1448 // track1->ApproximateHelixWithLine(mindist,GetFieldkG(),line1);
50ff0bcd 1449
1450 for(Int_t j=i+1; j<nacc; j++){
6d8df534 1451 track2 = (AliExternalTrackParam*)fTrkArraySel.At(j);
50ff0bcd 1452 alpha=track2->GetAlpha();
1453 mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
f09c879d 1454 track2->GetXYZAt(mindist,GetFieldkG(),pos);
1455 track2->GetPxPyPzAt(mindist,GetFieldkG(),dir);
50ff0bcd 1456 AliStrLine *line2 = new AliStrLine(pos,dir);
1457 // AliStrLine *line2 = new AliStrLine();
f09c879d 1458 // track2->ApproximateHelixWithLine(mindist,GetFieldkG(),line2);
50ff0bcd 1459 Double_t distCA=line2->GetDCA(line1);
fee9ef0d 1460 //printf("%d %d %f\n",i,j,distCA);
1461 if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
50ff0bcd 1462 Double_t pnt1[3],pnt2[3],crosspoint[3];
2d57349e 1463
50ff0bcd 1464 if(optUseWeights<=0){
1465 Int_t retcode = line2->Cross(line1,crosspoint);
1466 if(retcode>=0){
1467 ncombi++;
1468 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1469 for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1470 }
1471 }
1472 if(optUseWeights>0){
1473 Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2);
1474 if(retcode>=0){
5c03260f 1475 Double_t cs, sn;
50ff0bcd 1476 alpha=track1->GetAlpha();
1477 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1478 Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
1479 alpha=track2->GetAlpha();
1480 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1481 Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
1482 Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
1483 Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
1484 Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
1485 Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
1486 crosspoint[0]=wx1*pnt1[0] + wx2*pnt2[0];
1487 crosspoint[1]=wy1*pnt1[1] + wy2*pnt2[1];
1488 crosspoint[2]=wz1*pnt1[2] + wz2*pnt2[2];
1489
1490 ncombi++;
1491 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1492 for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1493 }
1494 }
1495 }
1496 delete line2;
1497 }
1498 delete line1;
71d84967 1499 }
50ff0bcd 1500 if(ncombi>0){
1501 for(Int_t jj=0;jj<3;jj++){
1502 initPos[jj] = aver[jj]/ncombi;
fee9ef0d 1503 //printf("%f\n",initPos[jj]);
50ff0bcd 1504 aversq[jj]/=ncombi;
1505 sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
1506 sigma+=sigmasq[jj];
1507 }
1508 sigma=TMath::Sqrt(TMath::Abs(sigma));
1509 }
1510 else {
1511 Warning("VertexFinder","Finder did not succed");
1512 sigma=999;
1513 }
1514 fVert.SetXYZ(initPos);
1515 fVert.SetDispersion(sigma);
1516 fVert.SetNContributors(ncombi);
2d57349e 1517}
ec1be5d5 1518//---------------------------------------------------------------------------
3f2db92f 1519void AliVertexerTracks::VertexFitter(Bool_t vfit, Bool_t chiCalc,Int_t useWeights)
fee9ef0d 1520{
ec1be5d5 1521//
1522// The optimal estimate of the vertex position is given by a "weighted
0bf2e2b4 1523// average of tracks positions".
fee9ef0d 1524// Original method: V. Karimaki, CMS Note 97/0051
ec1be5d5 1525//
3f2db92f 1526 const double kTiny = 1e-9;
6d8df534 1527 Bool_t useConstraint = fConstraint;
ec1be5d5 1528 Double_t initPos[3];
0bf2e2b4 1529 if(!fOnlyFitter) {
1530 fVert.GetXYZ(initPos);
1531 } else {
1532 initPos[0]=fNominalPos[0];
1533 initPos[1]=fNominalPos[1];
1534 initPos[2]=fNominalPos[2];
1535 }
1536
6d8df534 1537 Int_t nTrksSel = (Int_t)fTrkArraySel.GetEntries();
1538 if(nTrksSel==1) useConstraint=kTRUE;
3f2db92f 1539 AliDebug(1,Form("--- VertexFitter(): start (%d,%d,%d)",vfit,chiCalc,useWeights));
919e537f 1540 AliDebug(1,Form(" Number of tracks in array: %d\n",nTrksSel));
1541 AliDebug(1,Form(" Minimum # tracks required in fit: %d\n",fMinTracks));
1542 AliDebug(1,Form(" Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]));
1543 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 1544
0bf2e2b4 1545 // special treatment for few-tracks fits (e.g. secondary vertices)
3f2db92f 1546 Bool_t uUi3by3 = kFALSE; if(nTrksSel<5 && !useConstraint && !useWeights) uUi3by3 = kTRUE;
0bf2e2b4 1547
ec1be5d5 1548 Int_t i,j,k,step=0;
3f2db92f 1549 static TMatrixD rv(3,1);
1550 static TMatrixD vV(3,3);
ec1be5d5 1551 rv(0,0) = initPos[0];
1552 rv(1,0) = initPos[1];
3f2db92f 1553 rv(2,0) = initPos[2];
ec1be5d5 1554 Double_t xlStart,alpha;
3f2db92f 1555 Int_t nTrksUsed = 0;
1556 Double_t chi2=0,chi2i,chi2b;
6d8df534 1557 AliExternalTrackParam *t = 0;
ec1be5d5 1558 Int_t failed = 0;
1559
50ff0bcd 1560 // initial vertex covariance matrix
1561 TMatrixD vVb(3,3);
1562 vVb(0,0) = fNominalCov[0];
1563 vVb(0,1) = fNominalCov[1];
1564 vVb(0,2) = 0.;
1565 vVb(1,0) = fNominalCov[1];
1566 vVb(1,1) = fNominalCov[2];
1567 vVb(1,2) = 0.;
1568 vVb(2,0) = 0.;
1569 vVb(2,1) = 0.;
1570 vVb(2,2) = fNominalCov[5];
1571 TMatrixD vVbInv(TMatrixD::kInverted,vVb);
1572 TMatrixD rb(3,1);
1573 rb(0,0) = fNominalPos[0];
1574 rb(1,0) = fNominalPos[1];
1575 rb(2,0) = fNominalPos[2];
1576 TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb);
3f2db92f 1577 //
1578 int currBC = fVert.GetBC();
1579 //
eae1677e 1580 // 2 steps:
1581 // 1st - estimate of vtx using all tracks
1582 // 2nd - estimate of global chi2
3f2db92f 1583 //
eae1677e 1584 for(step=0; step<2; step++) {
3f2db92f 1585 if (step==0 && !vfit) continue;
1586 else if (step==1 && !chiCalc) continue;
ec1be5d5 1587 chi2 = 0.;
6d8df534 1588 nTrksUsed = 0;
3f2db92f 1589 fMVWSum = fMVWE2 = 0;
1590 if(step==1) { initPos[0]=rv(0,0); initPos[1]=rv(1,0); initPos[2]=rv(2,0);}
1591 AliDebug(2,Form("Step%d: inipos: %+f %+f %+f MinTr: %d, Sig2:%.2f)",step,initPos[0],initPos[1],initPos[2],fMinTracks,fMVSigma2));
1592 //
ec1be5d5 1593 TMatrixD sumWiri(3,1);
1594 TMatrixD sumWi(3,3);
1595 for(i=0; i<3; i++) {
1596 sumWiri(i,0) = 0.;
1597 for(j=0; j<3; j++) sumWi(j,i) = 0.;
1598 }
1599
fee9ef0d 1600 // mean vertex constraint
1601 if(useConstraint) {
50ff0bcd 1602 for(i=0;i<3;i++) {
1603 sumWiri(i,0) += vVbInvrb(i,0);
1604 for(k=0;k<3;k++) sumWi(i,k) += vVbInv(i,k);
07680cae 1605 }
fee9ef0d 1606 // chi2
1607 TMatrixD deltar = rv; deltar -= rb;
1608 TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar);
1609 chi2b = deltar(0,0)*vVbInvdeltar(0,0)+
1610 deltar(1,0)*vVbInvdeltar(1,0)+
1611 deltar(2,0)*vVbInvdeltar(2,0);
1612 chi2 += chi2b;
07680cae 1613 }
1614
ec1be5d5 1615 // loop on tracks
6d8df534 1616 for(k=0; k<nTrksSel; k++) {
3f2db92f 1617 //
ec1be5d5 1618 // get track from track array
6d8df534 1619 t = (AliExternalTrackParam*)fTrkArraySel.At(k);
3f2db92f 1620 if (useWeights && t->TestBit(kBitUsed)) continue;
1621 //
1622 int tBC = int(t->GetUniqueID()) - kTOFBCShift; // BC assigned to this track
1623 if (fSelectOnTOFBunchCrossing) {
1624 if (!fKeepAlsoUnflaggedTOFBunchCrossing) continue; // don't consider tracks with undefined BC
1625 if (currBC!=AliVTrack::kTOFBCNA && tBC!=AliVTrack::kTOFBCNA && tBC!=currBC) continue; // track does not match to current BCid
1626 }
ec1be5d5 1627 alpha = t->GetAlpha();
1628 xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
fee9ef0d 1629 // to vtxSeed (from finder)
6d8df534 1630 t->PropagateTo(xlStart,GetFieldkG());
fee9ef0d 1631
ec1be5d5 1632 // vector of track global coordinates
1633 TMatrixD ri(3,1);
fee9ef0d 1634 // covariance matrix of ri
1635 TMatrixD wWi(3,3);
1636
1637 // get space point from track
0bf2e2b4 1638 if(!TrackToPoint(t,ri,wWi,uUi3by3)) continue;
ec1be5d5 1639
1640 // track chi2
1641 TMatrixD deltar = rv; deltar -= ri;
1642 TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
1643 chi2i = deltar(0,0)*wWideltar(0,0)+
1644 deltar(1,0)*wWideltar(1,0)+
1645 deltar(2,0)*wWideltar(2,0);
3f2db92f 1646 //
1647 if (useWeights) {
1648 //double sg = TMath::Sqrt(fMVSigma2);
1649 //double chi2iw = (deltar(0,0)*wWideltar(0,0)+deltar(1,0)*wWideltar(1,0))/sg + deltar(2,0)*wWideltar(2,0)/fMVSigma2;
1650 //double wgh = (1-chi2iw/fMVTukey2);
1651 double chi2iw = chi2i;
1652 double wgh = (1-chi2iw/fMVTukey2/fMVSigma2);
1653
1654 if (wgh<kTiny) wgh = 0;
1655 else if (useWeights==2) wgh = 1.; // use as binary weight
1656 if (step==1) ((AliESDtrack*)t)->SetBit(kBitUsed, wgh>0);
1657 if (wgh<kTiny) continue; // discard the track
1658 wWi *= wgh; // RS: use weight?
1659 fMVWSum += wgh;
1660 fMVWE2 += wgh*chi2iw;
1661 }
ec1be5d5 1662 // add to total chi2
3f2db92f 1663 if (fSelectOnTOFBunchCrossing && tBC!=AliVTrack::kTOFBCNA && currBC<0) currBC = tBC;
1664 //
ec1be5d5 1665 chi2 += chi2i;
3f2db92f 1666 TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
ec1be5d5 1667 sumWiri += wWiri;
1668 sumWi += wWi;
1669
6d8df534 1670 nTrksUsed++;
ec1be5d5 1671 } // end loop on tracks
1672
6d8df534 1673 if(nTrksUsed < fMinTracks) {
ec1be5d5 1674 failed=1;
1675 continue;
1676 }
1677
1678 Double_t determinant = sumWi.Determinant();
f09c879d 1679 if(determinant < fMinDetFitter) {
919e537f 1680 AliDebug(1,Form("det(V) = %f (<%f)\n",determinant,fMinDetFitter));
ec1be5d5 1681 failed=1;
1682 continue;
1683 }
1684
0bf2e2b4 1685 if(step==0) {
1686 // inverted of weights matrix
1687 TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1688 vV = invsumWi;
1689 // position of primary vertex
1690 rv.Mult(vV,sumWiri);
1691 }
eae1677e 1692 } // end loop on the 2 steps
ec1be5d5 1693
ec1be5d5 1694 if(failed) {
3f2db92f 1695 fVert.SetNContributors(-1);
1696 if (chiCalc) {
1697 TooFewTracks();
1698 if (fCurrentVertex) fVert = *fCurrentVertex; // RS
1699 }
ec1be5d5 1700 return;
1701 }
3f2db92f 1702 //
ec1be5d5 1703 Double_t position[3];
1704 position[0] = rv(0,0);
1705 position[1] = rv(1,0);
1706 position[2] = rv(2,0);
1707 Double_t covmatrix[6];
1708 covmatrix[0] = vV(0,0);
1709 covmatrix[1] = vV(0,1);
1710 covmatrix[2] = vV(1,1);
1711 covmatrix[3] = vV(0,2);
1712 covmatrix[4] = vV(1,2);
1713 covmatrix[5] = vV(2,2);
1714
dc3719f3 1715 // for correct chi2/ndf, count constraint as additional "track"
6d8df534 1716 if(fConstraint) nTrksUsed++;
3f2db92f 1717 //
1718 if (vfit && !chiCalc) { // RS: special mode for multi-vertex finder
1719 fVert.SetXYZ(position);
1720 fVert.SetCovarianceMatrix(covmatrix);
1721 fVert.SetNContributors(nTrksUsed);
1722 return;
1723 }
ec1be5d5 1724 // store data in the vertex object
5b78b791 1725 if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
6d8df534 1726 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
3f2db92f 1727 fCurrentVertex->SetBC(currBC);
1728 fVert = *fCurrentVertex; // RS
919e537f 1729 AliDebug(1," Vertex after fit:");
1730 AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
1731 AliDebug(1,"--- VertexFitter(): finish\n");
1732
ec1be5d5 1733
1734 return;
1735}
50ff0bcd 1736//----------------------------------------------------------------------------
f5740e9a 1737AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(const TObjArray *trkArray,
6d8df534 1738 UShort_t *id,
1739 Bool_t optUseFitter,
97669588 1740 Bool_t optPropagate,
1741 Bool_t optUseDiamondConstraint)
fee9ef0d 1742{
eae1677e 1743//
6d8df534 1744// Return vertex from tracks (AliExternalTrackParam) in array
eae1677e 1745//
5b78b791 1746 fCurrentVertex = 0;
97669588 1747 // set optUseDiamondConstraint=TRUE only if you are reconstructing the
1748 // primary vertex!
1749 if(optUseDiamondConstraint) {
1750 SetConstraintOn();
1751 } else {
1752 SetConstraintOff();
1753 }
fee9ef0d 1754
50ff0bcd 1755 // get tracks and propagate them to initial vertex position
6d8df534 1756 fIdSel = new UShort_t[(Int_t)trkArray->GetEntriesFast()];
1757 Int_t nTrksSel = PrepareTracks(*trkArray,id,0);
97669588 1758 if((!optUseDiamondConstraint && nTrksSel<TMath::Max(2,fMinTracks)) ||
1759 (optUseDiamondConstraint && nTrksSel<1)) {
6d8df534 1760 TooFewTracks();
fee9ef0d 1761 return fCurrentVertex;
50ff0bcd 1762 }
1763
6d8df534 1764 // vertex finder
97669588 1765 if(nTrksSel==1) {
1766 AliDebug(1,"Just one track");
1767 OneTrackVertFinder();
1768 } else {
1769 switch (fAlgo) {
1770 case 1: StrLinVertexFinderMinDist(1); break;
1771 case 2: StrLinVertexFinderMinDist(0); break;
1772 case 3: HelixVertexFinder(); break;
1773 case 4: VertexFinder(1); break;
1774 case 5: VertexFinder(0); break;
1775 default: printf("Wrong algorithm\n"); break;
1776 }
fee9ef0d 1777 }
919e537f 1778 AliDebug(1," Vertex finding completed\n");
fee9ef0d 1779
1780 // vertex fitter
6d8df534 1781 if(optUseFitter) {
1782 VertexFitter();
1783 } else {
fee9ef0d 1784 Double_t position[3]={fVert.GetXv(),fVert.GetYv(),fVert.GetZv()};
146c29df 1785 Double_t covmatrix[6];
1786 fVert.GetCovMatrix(covmatrix);
fee9ef0d 1787 Double_t chi2=99999.;
6d8df534 1788 Int_t nTrksUsed=fVert.GetNContributors();
1789 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
fee9ef0d 1790 }
1791 fCurrentVertex->SetDispersion(fVert.GetDispersion());
1792
1793
1794 // set indices of used tracks and propagate track to found vertex
50ff0bcd 1795 UShort_t *indices = 0;
6d8df534 1796 Double_t d0z0[2],covd0z0[3];
1797 AliExternalTrackParam *t = 0;
fee9ef0d 1798 if(fCurrentVertex->GetNContributors()>0) {
44a76311 1799 indices = new UShort_t[fTrkArraySel.GetEntriesFast()];
6d8df534 1800 for(Int_t jj=0; jj<(Int_t)fTrkArraySel.GetEntriesFast(); jj++) {
1801 indices[jj] = fIdSel[jj];
1802 t = (AliExternalTrackParam*)fTrkArraySel.At(jj);
1803 if(optPropagate && optUseFitter) {
fee9ef0d 1804 if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
6d8df534 1805 t->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
919e537f 1806 AliDebug(1,Form("Track %d propagated to found vertex",jj));
6d8df534 1807 } else {
fee9ef0d 1808 AliWarning("Found vertex outside beam pipe!");
1809 }
1810 }
50ff0bcd 1811 }
fee9ef0d 1812 fCurrentVertex->SetIndices(fCurrentVertex->GetNContributors(),indices);
50ff0bcd 1813 }
6d8df534 1814
1815 // clean up
4ce766eb 1816 if (indices) {delete [] indices; indices=NULL;}
6fefa5ce 1817 delete [] fIdSel; fIdSel=NULL;
6d8df534 1818 fTrkArraySel.Delete();
50ff0bcd 1819
fee9ef0d 1820 return fCurrentVertex;
eae1677e 1821}
3f2db92f 1822
50ff0bcd 1823//----------------------------------------------------------------------------
6d8df534 1824AliESDVertex* AliVertexerTracks::VertexForSelectedESDTracks(TObjArray *trkArray,
3f2db92f 1825 Bool_t optUseFitter,
97669588 1826 Bool_t optPropagate,
1827 Bool_t optUseDiamondConstraint)
1828
fee9ef0d 1829{
50ff0bcd 1830//
6d8df534 1831// Return vertex from array of ESD tracks
50ff0bcd 1832//
5834e9a1 1833
6d8df534 1834 Int_t nTrks = (Int_t)trkArray->GetEntriesFast();
1835 UShort_t *id = new UShort_t[nTrks];
1836
1837 AliESDtrack *esdt = 0;
50ff0bcd 1838 for(Int_t i=0; i<nTrks; i++){
6d8df534 1839 esdt = (AliESDtrack*)trkArray->At(i);
1840 id[i] = (UShort_t)esdt->GetID();
50ff0bcd 1841 }
1842
97669588 1843 VertexForSelectedTracks(trkArray,id,optUseFitter,optPropagate,optUseDiamondConstraint);
6d8df534 1844
4ce766eb 1845 delete [] id; id=NULL;
6d8df534 1846
1847 return fCurrentVertex;
50ff0bcd 1848}
3f2db92f 1849
1850//______________________________________________________
1851Bool_t AliVertexerTracks::FindNextVertexMV()
1852{
1853 // try to find a new vertex
1854 fMVSigma2 = fMVSig2Ini;
1855 double prevSig2 = fMVSigma2;
1856 double minDst2 = fMVMinDst*fMVMinDst;
1857 const double kSigLimit = 1.0;
1858 const double kSigLimitE = kSigLimit+1e-6;
1859 const double kPushFactor = 0.5;
1860 const int kMaxIter = 20;
1861 double push = kPushFactor;
1862 //
1863 int iter = 0;
1864 double posP[3]={0,0,0},pos[3]={0,0,0};
1865 fVert.GetXYZ(posP);
1866 //
1867
1868 do {
1869 fVert.SetBC(AliVTrack::kTOFBCNA);
1870 VertexFitter(kTRUE,kFALSE,1);
1871 if (fVert.GetNContributors()<fMinTracks) {
1872 AliDebug(3,Form("Failed in iteration %d: No Contributirs",iter));
1873 break;
1874 } // failed
1875 if (fMVWSum>0) fMVSigma2 = TMath::Max(fMVWE2/fMVWSum,kSigLimit);
1876 else {
1877 AliDebug(3,Form("Failed %d, no weithgs",iter));
1878 iter = kMaxIter+1;
1879 break;
1880 } // failed
1881 //
1882 double sigRed = (prevSig2-fMVSigma2)/prevSig2; // reduction of sigma2
1883 //
1884 fVert.GetXYZ(pos);
1885 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]);
1886 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));
1887 if ( (++iter<kMaxIter) && (sigRed<0 || sigRed<fMVMinSig2Red) && fMVSigma2>fMVMaxSigma2) {
1888 fMVSigma2 *= push; // stuck, push little bit
1889 push *= kPushFactor;
1890 if (fMVSigma2<1.) fMVSigma2 = 1.;
1891 AliDebug(3,Form("Pushed sigma2 to %f",fMVSigma2));
1892 }
1893 else if (dst2<minDst2 && ((sigRed<0 || sigRed<fMVMinSig2Red) || fMVSigma2<kSigLimitE)) break;
1894 //
1895 fVert.GetXYZ(posP); // fetch previous vertex position
1896 prevSig2 = fMVSigma2;
1897 } while(iter<kMaxIter);
1898 //
1899 if (fVert.GetNContributors()<0 || iter>kMaxIter || fMVSigma2>fMVMaxSigma2) {
1900 return kFALSE;
1901 }
1902 else {
1903 VertexFitter(kFALSE,kTRUE,fMVFinalWBinary ? 2:1); // final chi2 calculation
1904 int nv = fMVVertices->GetEntries();
1905 // create indices
1906 int ntrk = fTrkArraySel.GetEntries();
1907 int nindices = fCurrentVertex->GetNContributors() - (fConstraint ? 1:0);
1908 UShort_t *indices = 0;
1909 if (nindices>0) indices = new UShort_t[nindices];
1910 int nadded = 0;
1911 for (int itr=0;itr<ntrk;itr++) {
1912 AliExternalTrackParam* t = (AliExternalTrackParam*)fTrkArraySel[itr];
1913 if (t->TestBit(kBitAccounted) || !t->TestBit(kBitUsed)) continue; // already belongs to some vertex
1914 t->SetBit(kBitAccounted);
1915 indices[nadded++] = fIdSel[itr];
1916 }
1917 if (nadded!=nindices) printf("Mismatch : NInd: %d Nadd: %d\n",nindices,nadded);
1918 fCurrentVertex->SetIndices(nindices,indices);
1919 // set vertex title
bac16cc1 1920 TString title="VertexerTracksMVNoConstraint";
1921 if(fConstraint) title="VertexerTracksMVWithConstraint";
3f2db92f 1922 fCurrentVertex->SetTitle(title.Data());
1923 fMVVertices->AddLast(fCurrentVertex);
1924 AliDebug(3,Form("Added new vertex #%d NCont:%d XYZ: %f %f %f",nindices,nv,fCurrentVertex->GetX(),fCurrentVertex->GetY(),fCurrentVertex->GetZ()));
1925 if (indices) delete[] indices;
1926 fCurrentVertex = 0; // already attached to fMVVertices
1927 return kTRUE;
1928 }
1929}
1930
1931//______________________________________________________
1932void AliVertexerTracks::FindVerticesMV()
1933{
1934 // find and fit multiple vertices
1935 //
1936 double step = fMVScanStep>1 ? fMVScanStep : 1.;
1937 double zmx = 3*TMath::Sqrt(fNominalCov[5]);
1938 double zmn = -zmx;
1939 int nz = (zmx-zmn)/step; if (nz<1) nz=1;
1940 double dz = (zmx-zmn)/nz;
1941 int izStart=0;
1942 AliDebug(2,Form("%d seeds between %f and %f",nz,zmn+dz/2,zmx+dz/2));
1943 //
1944 if (!fMVVertices) fMVVertices = new TObjArray(10);
1945 fMVVertices->Clear();
1946 //
1947 int ntrLeft = (Int_t)fTrkArraySel.GetEntries();
1948 //
1949 double sig2Scan = fMVSig2Ini;
1950 Bool_t runMore = kTRUE;
1951 int cntWide = 0;
1952 while (runMore) {
1953 fMVSig2Ini = sig2Scan*1e3; // try wide search
1954 Bool_t found = kFALSE;
1955 cntWide++;
1956 fVert.SetNContributors(-1);
1957 fVert.SetXYZ(fNominalPos);
1958 AliDebug(3,Form("Wide search #%d Z= %f Sigma2=%f",cntWide,fNominalPos[2],fMVSig2Ini));
1959 if (FindNextVertexMV()) { // are there tracks left to consider?
1960 AliESDVertex* vtLast = (AliESDVertex*)fMVVertices->Last();
1961 if (vtLast && vtLast->GetNContributors()>0) ntrLeft -= vtLast->GetNContributors()-(fConstraint ? 1:0);
1962 if (ntrLeft<1) runMore = kFALSE;
1963 found = kTRUE;
1964 continue;
1965 }
1966 // if nothing is found, do narrow sig2ini scan
1967 fMVSig2Ini = sig2Scan;
1968 for (;izStart<nz;izStart++) {
1969 double zSeed = zmn+dz*(izStart+0.5);
1970 AliDebug(3,Form("Seed %d: Z= %f Sigma2=%f",izStart,zSeed,fMVSig2Ini));
1971 fVert.SetNContributors(-1);
1972 fVert.SetXYZ(fNominalPos);
1973 fVert.SetZv(zSeed);
1974 if (FindNextVertexMV()) { // are there tracks left to consider?
1975 AliESDVertex* vtLast = (AliESDVertex*)fMVVertices->Last();
1976 if (vtLast && vtLast->GetNContributors()>0) ntrLeft -= vtLast->GetNContributors()-(fConstraint ? 1:0);
1977 if (ntrLeft<1) runMore = kFALSE;
1978 found = kTRUE;
1979 break;
1980 }
1981 }
1982 runMore = found; // if nothing was found, no need for new iteration
1983 }
1984 fMVSig2Ini = sig2Scan;
1985 int nvFound = fMVVertices->GetEntriesFast();
1986 AliDebug(2,Form("Number of found vertices: %d",nvFound));
1987 if (nvFound<1) TooFewTracks();
1988 if (AliLog::GetGlobalDebugLevel()>0) fMVVertices->Print();
1989 //
1990}
1991
1992//______________________________________________________
1993void AliVertexerTracks::AnalyzePileUp(AliESDEvent* esdEv)
1994{
1995 // if multiple vertices are found, try to find the primary one and attach it as fCurrentVertex
1996 // then attach pile-up vertices directly to esdEv
1997 //
1998 int nFND = (fMVVertices && fMVVertices->GetEntriesFast()) ? fMVVertices->GetEntriesFast() : 0;
1999 if (nFND<1) { if (!fCurrentVertex) TooFewTracks(); return;} // no multiple vertices
2000 //
2001 int indCont[nFND];
2002 int nIndx[nFND];
2003 for (int iv=0;iv<nFND;iv++) {
2004 AliESDVertex* fnd = (AliESDVertex*)fMVVertices->At(iv);
2005 indCont[iv] = iv;
2006 nIndx[iv] = fnd->GetNIndices();
2007 }
2008 TMath::Sort(nFND, nIndx, indCont, kTRUE); // sort in decreasing order of Nindices
2009 double dists[nFND];
2010 int distOrd[nFND],indx[nFND];
2011 for (int iv=0;iv<nFND;iv++) {
2012 AliESDVertex* fndI = (AliESDVertex*)fMVVertices->At(indCont[iv]);
2013 if (fndI->GetStatus()<1) continue; // discarded
2014 int ncomp = 0;
2015 for (int jv=iv+1;jv<nFND;jv++) {
2016 AliESDVertex* fndJ = (AliESDVertex*)fMVVertices->At(indCont[jv]);
2017 if (fndJ->GetStatus()<1) continue;
2018 dists[ncomp] = fndI->GetWDist(fndJ)*fndJ->GetNIndices();
2019 distOrd[ncomp] = indCont[jv];
2020 indx[ncomp] = ncomp;
2021 ncomp++;
2022 }
2023 if (ncomp<1) continue;
2024 TMath::Sort(ncomp, dists, indx, kFALSE); // sort in increasing distance order
2025 for (int jv0=0;jv0<ncomp;jv0++) {
2026 int jv = distOrd[indx[jv0]];
2027 AliESDVertex* fndJ = (AliESDVertex*)fMVVertices->At(jv);
2028 if (dists[indx[jv0]]<fMVMaxWghNtr) { // candidate for split vertex
2029 //before eliminating the small close vertex, check if we should transfere its BC to largest one
2030 if (fndJ->GetBC()!=AliVTrack::kTOFBCNA && fndI->GetBC()==AliVTrack::kTOFBCNA) fndI->SetBC(fndJ->GetBC());
2031 //
2032 // leave the vertex only if both vertices have definit but different BC's, then this is not splitting.
2033 if ( (fndJ->GetBC()==fndI->GetBC()) || (fndJ->GetBC()==AliVTrack::kTOFBCNA)) fndJ->SetNContributors(-fndJ->GetNContributors());
2034 }
2035 }
2036 }
2037 //
2038 // select as a primary the largest vertex with BC=0, or the largest with BC non-ID
2039 int primBC0=-1,primNoBC=-1;
2040 for (int iv0=0;iv0<nFND;iv0++) {
2041 int iv = indCont[iv0];
2042 AliESDVertex* fndI = (AliESDVertex*)fMVVertices->At(iv);
2043 if (!fndI) continue;
2044 if (fndI->GetStatus()<1) {fMVVertices->RemoveAt(iv); delete fndI; continue;}
2045 if (primBC0<0 && fndI->GetBC()==0) primBC0 = iv;
2046 if (primNoBC<0 && fndI->GetBC()==AliVTrack::kTOFBCNA) primNoBC = iv;
2047 }
2048 //
2049 if (primBC0>=0) fCurrentVertex = (AliESDVertex*)fMVVertices->At(primBC0);
2050 if (!fCurrentVertex && primNoBC>=0) fCurrentVertex = (AliESDVertex*)fMVVertices->At(primNoBC);
2051 if (fCurrentVertex) fMVVertices->Remove(fCurrentVertex);
2052 else { // all vertices have BC>0, no primary vertex
2053 fCurrentVertex = new AliESDVertex();
2054 fCurrentVertex->SetNContributors(-3);
2055 fCurrentVertex->SetBC(AliVTrack::kTOFBCNA);
2056 }
2057 fCurrentVertex->SetID(-1);
2058 //
2059 // add pileup vertices
2060 int nadd = 0;
2061 for (int iv0=0;iv0<nFND;iv0++) {
2062 int iv = indCont[iv0];
2063 AliESDVertex* fndI = (AliESDVertex*)fMVVertices->At(iv);
2064 if (!fndI) continue;
2065 fndI->SetID(++nadd);
2066 esdEv->AddPileupVertexTracks(fndI);
2067 }
2068 //
2069 fMVVertices->Delete();
2070 //
2071}
2072