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