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