Multi-vertexerTracks with clusterization + steering options in GRPRecoParam
[u/mrichter/AliRoot.git] / STEER / ESD / AliVertexerTracks.cxx
CommitLineData
2d57349e 1/**************************************************************************
2 * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//-----------------------------------------------------------------
6d8df534 17// Implementation of the vertexer from tracks
2d57349e 18//
19// Origin: AliITSVertexerTracks
20// A.Dainese, Padova,
21// andrea.dainese@pd.infn.it
22// M.Masera, Torino,
23// massimo.masera@to.infn.it
24// Moved to STEER and adapted to ESD tracks:
6d8df534 25// F.Prino, Torino, prino@to.infn.it
2d57349e 26//-----------------------------------------------------------------
27
28//---- Root headers --------
817e1004 29#include <TSystem.h>
30#include <TClonesArray.h>
a7ede274 31#include <TDirectory.h>
32#include <TFile.h>
2d57349e 33//---- AliRoot headers -----
34#include "AliStrLine.h"
6d8df534 35#include "AliExternalTrackParam.h"
2258e165 36#include "AliNeutralTrackParam.h"
0aa2faf4 37#include "AliVEvent.h"
38#include "AliVTrack.h"
2d57349e 39#include "AliESDtrack.h"
3f2db92f 40#include "AliESDEvent.h"
6d8df534 41#include "AliVertexerTracks.h"
2d57349e 42
3f2db92f 43
2d57349e 44ClassImp(AliVertexerTracks)
45
46
47//----------------------------------------------------------------------------
48AliVertexerTracks::AliVertexerTracks():
07680cae 49TObject(),
50fVert(),
51fCurrentVertex(0),
f09c879d 52fMode(0),
dc3719f3 53fFieldkG(-999.),
6d8df534 54fTrkArraySel(),
55fIdSel(0),
07680cae 56fTrksToSkip(0),
57fNTrksToSkip(0),
f09c879d 58fConstraint(kFALSE),
59fOnlyFitter(kFALSE),
60fMinTracks(1),
6b29d399 61fMinClusters(5),
f09c879d 62fDCAcut(0.1),
63fDCAcutIter0(0.1),
64fNSigma(3.),
fee9ef0d 65fMaxd0z0(0.5),
f09c879d 66fMinDetFitter(100.),
67fMaxTgl(1000.),
50ff0bcd 68fITSrefit(kTRUE),
40cfa14d 69fITSpureSA(kFALSE),
6b29d399 70fFiducialR(3.),
71fFiducialZ(30.),
3b768a1d 72fnSigmaForUi00(1.5),
8c75f668 73fAlgo(1),
ce1c94a7 74fAlgoIter0(4),
75fSelectOnTOFBunchCrossing(kFALSE),
3f2db92f 76fKeepAlsoUnflaggedTOFBunchCrossing(kTRUE),
77fMVWSum(0),
78fMVWE2(0),
79fMVTukey2(6.),
80fMVSigma2(1.),
81fMVSig2Ini(1.e3),
82fMVMaxSigma2(3.),
83fMVMinSig2Red(0.005),
84fMVMinDst(10.e-4),
85fMVScanStep(3.),
86fMVMaxWghNtr(10.),
87fMVFinalWBinary(kTRUE),
88fBCSpacing(50),
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),
6b29d399 114fMinClusters(5),
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;
333 default: printf("Wrong algorithm\n"); break;
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;
919e537f 387 default: printf("Wrong algorithm"); 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++) {
6d8df534 884 if(inindices[k]==id[l]) copyindex=kFALSE;
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();
1014 if (ncuts>4) SetMinClusters((Int_t)(TMath::Abs(cuts[3])));
1015 if (ncuts>5) SetMinTracks((Int_t)(cuts[4]));
1016 if (ncuts>6) SetNSigmad0(cuts[5]);
1017 if (ncuts>7) SetMinDetFitter(cuts[6]);
1018 if (ncuts>8) SetMaxTgl(cuts[7]);
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)
dc3719f3 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);
8c75f668 1152 if(!TrackToPoint(track1,ri,wWi)) {optUseWeights=kFALSE;printf("WARNING\n");}
146c29df 1153 Double_t wmat[9];
1154 Int_t iel=0;
1155 for(Int_t ia=0;ia<3;ia++){
1156 for(Int_t ib=0;ib<3;ib++){
1157 wmat[iel]=wWi(ia,ib);
1158 iel++;
1159 }
1160 }
8e104736 1161 linarray[i] = new AliStrLine(pos,sigmasq,wmat,dir);
dc3719f3 1162 }
8e104736 1163 fVert=TrackletVertexFinder(linarray,knacc,optUseWeights);
1164 for(Int_t i=0; i<knacc; i++) delete linarray[i];
1165 delete [] linarray;
dc3719f3 1166}
1167//---------------------------------------------------------------------------
f5740e9a 1168AliESDVertex AliVertexerTracks::TrackletVertexFinder(const TClonesArray *lines, Int_t optUseWeights)
fee9ef0d 1169{
8e104736 1170 // Calculate the point at minimum distance to prepared tracks (TClonesArray)
dc3719f3 1171 const Int_t knacc = (Int_t)lines->GetEntriesFast();
8e104736 1172 AliStrLine** lines2 = new AliStrLine* [knacc];
1173 for(Int_t i=0; i<knacc; i++){
1174 lines2[i]= (AliStrLine*)lines->At(i);
1175 }
1176 AliESDVertex vert = TrackletVertexFinder(lines2,knacc,optUseWeights);
1177 delete [] lines2;
1178 return vert;
1179}
1180
1181//---------------------------------------------------------------------------
1182AliESDVertex AliVertexerTracks::TrackletVertexFinder(AliStrLine **lines, const Int_t knacc, Int_t optUseWeights)
1183{
1184 // Calculate the point at minimum distance to prepared tracks (array of AliStrLine)
1185
dc3719f3 1186 Double_t initPos[3]={0.,0.,0.};
2d57349e 1187
3f2db92f 1188 Double_t (*vectP0)[3]=new Double_t [knacc][3];
1189 Double_t (*vectP1)[3]=new Double_t [knacc][3];
2d57349e 1190
1191 Double_t sum[3][3];
1192 Double_t dsum[3]={0,0,0};
146c29df 1193 TMatrixD sumWi(3,3);
1194 for(Int_t i=0;i<3;i++){
1195 for(Int_t j=0;j<3;j++){
1196 sum[i][j]=0;
1197 sumWi(i,j)=0.;
1198 }
1199 }
1200
2d57349e 1201 for(Int_t i=0; i<knacc; i++){
8e104736 1202 AliStrLine *line1 = lines[i];
dc3719f3 1203 Double_t p0[3],cd[3],sigmasq[3];
146c29df 1204 Double_t wmat[9];
1100f3cd 1205 if(!line1) {printf("ERROR %d %d\n",i,knacc); continue;}
2d57349e 1206 line1->GetP0(p0);
1207 line1->GetCd(cd);
dc3719f3 1208 line1->GetSigma2P0(sigmasq);
146c29df 1209 line1->GetWMatrix(wmat);
1210 TMatrixD wWi(3,3);
1211 Int_t iel=0;
1212 for(Int_t ia=0;ia<3;ia++){
1213 for(Int_t ib=0;ib<3;ib++){
1214 wWi(ia,ib)=wmat[iel];
1215 iel++;
1216 }
1217 }
1218
1219 sumWi+=wWi;
1220
2d57349e 1221 Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
1222 vectP0[i][0]=p0[0];
1223 vectP0[i][1]=p0[1];
1224 vectP0[i][2]=p0[2];
1225 vectP1[i][0]=p1[0];
1226 vectP1[i][1]=p1[1];
1227 vectP1[i][2]=p1[2];
1228
1229 Double_t matr[3][3];
1230 Double_t dknow[3];
1231 if(optUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
dc3719f3 1232 else GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
1233
2d57349e 1234
1235 for(Int_t iii=0;iii<3;iii++){
1236 dsum[iii]+=dknow[iii];
1237 for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
1238 }
2d57349e 1239 }
146c29df 1240
1241 TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1242 Double_t covmatrix[6];
1243 covmatrix[0] = invsumWi(0,0);
1244 covmatrix[1] = invsumWi(0,1);
1245 covmatrix[2] = invsumWi(1,1);
1246 covmatrix[3] = invsumWi(0,2);
1247 covmatrix[4] = invsumWi(1,2);
1248 covmatrix[5] = invsumWi(2,2);
1249
2d57349e 1250 Double_t vett[3][3];
1251 Double_t det=GetDeterminant3X3(sum);
dc3719f3 1252 Double_t sigma=0;
2d57349e 1253
f5740e9a 1254 if(TMath::Abs(det) > kAlmost0){
fee9ef0d 1255 for(Int_t zz=0;zz<3;zz++){
1256 for(Int_t ww=0;ww<3;ww++){
1257 for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
1258 }
1259 for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
1260 initPos[zz]=GetDeterminant3X3(vett)/det;
1261 }
1262
1263
1264 for(Int_t i=0; i<knacc; i++){
1265 Double_t p0[3]={0,0,0},p1[3]={0,0,0};
1266 for(Int_t ii=0;ii<3;ii++){
1267 p0[ii]=vectP0[i][ii];
1268 p1[ii]=vectP1[i][ii];
1269 }
1270 sigma+=GetStrLinMinDist(p0,p1,initPos);
1271 }
1272
f09c879d 1273 if(sigma>0.) {sigma=TMath::Sqrt(sigma);}else{sigma=999;}
fee9ef0d 1274 }else{
50ff0bcd 1275 sigma=999;
1276 }
146c29df 1277 AliESDVertex theVert(initPos,covmatrix,99999.,knacc);
1278 theVert.SetDispersion(sigma);
4ce766eb 1279 delete [] vectP0;
1280 delete [] vectP1;
dc3719f3 1281 return theVert;
50ff0bcd 1282}
8e104736 1283
50ff0bcd 1284//---------------------------------------------------------------------------
6d8df534 1285Bool_t AliVertexerTracks::TrackToPoint(AliExternalTrackParam *t,
817e1004 1286 TMatrixD &ri,TMatrixD &wWi,
0bf2e2b4 1287 Bool_t uUi3by3) const
fee9ef0d 1288{
1289//
6d8df534 1290// Extract from the AliExternalTrackParam the global coordinates ri and covariance matrix
fee9ef0d 1291// wWi of the space point that it represents (to be used in VertexFitter())
1292//
1293
1294
1295 Double_t rotAngle = t->GetAlpha();
1296 if(rotAngle<0.) rotAngle += 2.*TMath::Pi();
1297 Double_t cosRot = TMath::Cos(rotAngle);
1298 Double_t sinRot = TMath::Sin(rotAngle);
3f2db92f 1299 /*
1300 // RS >>>
1301 Double_t lambda = TMath::ATan(t->GetTgl());
1302 Double_t cosLam = TMath::Cos(lambda);
1303 Double_t sinLam = TMath::Sin(lambda);
1304 // RS <<<
1305 */
fee9ef0d 1306 ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
1307 ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
1308 ri(2,0) = t->GetZ();
1309
0bf2e2b4 1310 if(!uUi3by3) {
1311 // matrix to go from global (x,y,z) to local (y,z);
1312 TMatrixD qQi(2,3);
1313 qQi(0,0) = -sinRot;
1314 qQi(0,1) = cosRot;
1315 qQi(0,2) = 0.;
3f2db92f 1316 //
0bf2e2b4 1317 qQi(1,0) = 0.;
1318 qQi(1,1) = 0.;
1319 qQi(1,2) = 1.;
3f2db92f 1320 //
1321 // RS: Added polar inclination
1322 /*
1323 qQi(1,0) = -sinLam*cosRot;
1324 qQi(1,1) = -sinLam*sinRot;
1325 qQi(1,2) = cosLam;
1326 */
0bf2e2b4 1327 // covariance matrix of local (y,z) - inverted
0bf2e2b4 1328 TMatrixD uUi(2,2);
6d8df534 1329 uUi(0,0) = t->GetSigmaY2();
1330 uUi(0,1) = t->GetSigmaZY();
1331 uUi(1,0) = t->GetSigmaZY();
1332 uUi(1,1) = t->GetSigmaZ2();
919e537f 1333 //printf(" Ui :");
1334 //printf(" %f %f",uUi(0,0),uUi(0,1));
1335 //printf(" %f %f",uUi(1,0),uUi(1,1));
0bf2e2b4 1336
1337 if(uUi.Determinant() <= 0.) return kFALSE;
1338 TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1339
1340 // weights matrix: wWi = qQiT * uUiInv * qQi
1341 TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1342 TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1343 wWi = m;
1344 } else {
1345 if(fVert.GetNContributors()<1) AliFatal("Vertex from finder is empty");
1346 // matrix to go from global (x,y,z) to local (x,y,z);
1347 TMatrixD qQi(3,3);
1348 qQi(0,0) = cosRot;
1349 qQi(0,1) = sinRot;
1350 qQi(0,2) = 0.;
1351 qQi(1,0) = -sinRot;
1352 qQi(1,1) = cosRot;
1353 qQi(1,2) = 0.;
1354 qQi(2,0) = 0.;
1355 qQi(2,1) = 0.;
1356 qQi(2,2) = 1.;
1357
1358 // covariance of fVert along the track
1359 Double_t p[3],pt,ptot;
1360 t->GetPxPyPz(p);
1361 pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
1362 ptot = TMath::Sqrt(pt*pt+p[2]*p[2]);
3b768a1d 1363 Double_t cphi = p[0]/pt; //cos(phi)=px/pt
1364 Double_t sphi = p[1]/pt; //sin(phi)=py/pt
1365 Double_t clambda = pt/ptot; //cos(lambda)=pt/ptot
1366 Double_t slambda = p[2]/ptot; //sin(lambda)=pz/ptot
0bf2e2b4 1367 Double_t covfVert[6];
1368 fVert.GetCovMatrix(covfVert);
1369 Double_t covfVertalongt =
3b768a1d 1370 covfVert[0]*cphi*cphi*clambda*clambda
1371 +covfVert[1]*2.*cphi*sphi*clambda*clambda
1372 +covfVert[3]*2.*cphi*clambda*slambda
1373 +covfVert[2]*sphi*sphi*clambda*clambda
1374 +covfVert[4]*2.*sphi*clambda*slambda
1375 +covfVert[5]*slambda*slambda;
0bf2e2b4 1376 // covariance matrix of local (x,y,z) - inverted
1377 TMatrixD uUi(3,3);
3b768a1d 1378 uUi(0,0) = covfVertalongt * fnSigmaForUi00 * fnSigmaForUi00;
919e537f 1379 AliDebug(1,Form("=====> sqrtUi00 cm %f",TMath::Sqrt(uUi(0,0))));
0bf2e2b4 1380 uUi(0,1) = 0.;
1381 uUi(0,2) = 0.;
1382 uUi(1,0) = 0.;
6d8df534 1383 uUi(1,1) = t->GetSigmaY2();
1384 uUi(1,2) = t->GetSigmaZY();
0bf2e2b4 1385 uUi(2,0) = 0.;
6d8df534 1386 uUi(2,1) = t->GetSigmaZY();
1387 uUi(2,2) = t->GetSigmaZ2();
0bf2e2b4 1388 //printf(" Ui :\n");
1389 //printf(" %f %f\n",uUi(0,0),uUi(0,1));
1390 //printf(" %f %f\n",uUi(1,0),uUi(1,1));
fee9ef0d 1391
0bf2e2b4 1392 if(uUi.Determinant() <= 0.) return kFALSE;
1393 TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1394
1395 // weights matrix: wWi = qQiT * uUiInv * qQi
1396 TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1397 TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1398 wWi = m;
1399 }
1400
fee9ef0d 1401
1402 return kTRUE;
1403}
1404//---------------------------------------------------------------------------
6d8df534 1405void AliVertexerTracks::TooFewTracks()
fee9ef0d 1406{
50ff0bcd 1407//
6d8df534 1408// When the number of tracks is < fMinTracks,
1409// deal with vertices not found and prepare to exit
50ff0bcd 1410//
919e537f 1411 AliDebug(1,"TooFewTracks");
2d57349e 1412
50ff0bcd 1413 Double_t pos[3],err[3];
50ff0bcd 1414 pos[0] = fNominalPos[0];
1415 err[0] = TMath::Sqrt(fNominalCov[0]);
1416 pos[1] = fNominalPos[1];
1417 err[1] = TMath::Sqrt(fNominalCov[2]);
6d8df534 1418 pos[2] = fNominalPos[2];
1419 err[2] = TMath::Sqrt(fNominalCov[5]);
1420 Int_t ncontr = (err[0]>1. ? -1 : -3);
5b78b791 1421 if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
50ff0bcd 1422 fCurrentVertex = new AliESDVertex(pos,err);
1423 fCurrentVertex->SetNContributors(ncontr);
2d57349e 1424
fee9ef0d 1425 if(fConstraint) {
1426 fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
1427 } else {
50ff0bcd 1428 fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
fee9ef0d 1429 }
2d57349e 1430
6d8df534 1431 if(!fTrkArraySel.IsEmpty()) fTrkArraySel.Delete();
6fefa5ce 1432 if(fIdSel) {delete [] fIdSel; fIdSel=NULL;}
1433 if(fTrksToSkip) {delete [] fTrksToSkip; fTrksToSkip=NULL;}
6d8df534 1434
50ff0bcd 1435 return;
1436}
1437//---------------------------------------------------------------------------
fee9ef0d 1438void AliVertexerTracks::VertexFinder(Int_t optUseWeights)
1439{
2d57349e 1440
50ff0bcd 1441 // Get estimate of vertex position in (x,y) from tracks DCA
1442
1443 Double_t initPos[3];
1444 initPos[2] = 0.;
1445 for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
6d8df534 1446 Int_t nacc = (Int_t)fTrkArraySel.GetEntriesFast();
50ff0bcd 1447 Double_t aver[3]={0.,0.,0.};
1448 Double_t aversq[3]={0.,0.,0.};
1449 Double_t sigmasq[3]={0.,0.,0.};
1450 Double_t sigma=0;
1451 Int_t ncombi = 0;
6d8df534 1452 AliExternalTrackParam *track1;
1453 AliExternalTrackParam *track2;
50ff0bcd 1454 Double_t pos[3],dir[3];
1455 Double_t alpha,mindist;
2d57349e 1456
50ff0bcd 1457 for(Int_t i=0; i<nacc; i++){
6d8df534 1458 track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
50ff0bcd 1459 alpha=track1->GetAlpha();
1460 mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
f09c879d 1461 track1->GetXYZAt(mindist,GetFieldkG(),pos);
1462 track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
50ff0bcd 1463 AliStrLine *line1 = new AliStrLine(pos,dir);
2d57349e 1464
50ff0bcd 1465 // AliStrLine *line1 = new AliStrLine();
f09c879d 1466 // track1->ApproximateHelixWithLine(mindist,GetFieldkG(),line1);
50ff0bcd 1467
1468 for(Int_t j=i+1; j<nacc; j++){
6d8df534 1469 track2 = (AliExternalTrackParam*)fTrkArraySel.At(j);
50ff0bcd 1470 alpha=track2->GetAlpha();
1471 mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
f09c879d 1472 track2->GetXYZAt(mindist,GetFieldkG(),pos);
1473 track2->GetPxPyPzAt(mindist,GetFieldkG(),dir);
50ff0bcd 1474 AliStrLine *line2 = new AliStrLine(pos,dir);
1475 // AliStrLine *line2 = new AliStrLine();
f09c879d 1476 // track2->ApproximateHelixWithLine(mindist,GetFieldkG(),line2);
50ff0bcd 1477 Double_t distCA=line2->GetDCA(line1);
fee9ef0d 1478 //printf("%d %d %f\n",i,j,distCA);
1479 if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
50ff0bcd 1480 Double_t pnt1[3],pnt2[3],crosspoint[3];
2d57349e 1481
50ff0bcd 1482 if(optUseWeights<=0){
1483 Int_t retcode = line2->Cross(line1,crosspoint);
1484 if(retcode>=0){
1485 ncombi++;
1486 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1487 for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1488 }
1489 }
1490 if(optUseWeights>0){
1491 Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2);
1492 if(retcode>=0){
5c03260f 1493 Double_t cs, sn;
50ff0bcd 1494 alpha=track1->GetAlpha();
1495 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1496 Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
1497 alpha=track2->GetAlpha();
1498 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1499 Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
1500 Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
1501 Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
1502 Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
1503 Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
1504 crosspoint[0]=wx1*pnt1[0] + wx2*pnt2[0];
1505 crosspoint[1]=wy1*pnt1[1] + wy2*pnt2[1];
1506 crosspoint[2]=wz1*pnt1[2] + wz2*pnt2[2];
1507
1508 ncombi++;
1509 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1510 for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1511 }
1512 }
1513 }
1514 delete line2;
1515 }
1516 delete line1;
71d84967 1517 }
50ff0bcd 1518 if(ncombi>0){
1519 for(Int_t jj=0;jj<3;jj++){
1520 initPos[jj] = aver[jj]/ncombi;
fee9ef0d 1521 //printf("%f\n",initPos[jj]);
50ff0bcd 1522 aversq[jj]/=ncombi;
1523 sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
1524 sigma+=sigmasq[jj];
1525 }
1526 sigma=TMath::Sqrt(TMath::Abs(sigma));
1527 }
1528 else {
1529 Warning("VertexFinder","Finder did not succed");
1530 sigma=999;
1531 }
1532 fVert.SetXYZ(initPos);
1533 fVert.SetDispersion(sigma);
1534 fVert.SetNContributors(ncombi);
2d57349e 1535}
ec1be5d5 1536//---------------------------------------------------------------------------
3f2db92f 1537void AliVertexerTracks::VertexFitter(Bool_t vfit, Bool_t chiCalc,Int_t useWeights)
fee9ef0d 1538{
ec1be5d5 1539//
1540// The optimal estimate of the vertex position is given by a "weighted
0bf2e2b4 1541// average of tracks positions".
fee9ef0d 1542// Original method: V. Karimaki, CMS Note 97/0051
ec1be5d5 1543//
3f2db92f 1544 const double kTiny = 1e-9;
6d8df534 1545 Bool_t useConstraint = fConstraint;
ec1be5d5 1546 Double_t initPos[3];
0bf2e2b4 1547 if(!fOnlyFitter) {
1548 fVert.GetXYZ(initPos);
1549 } else {
1550 initPos[0]=fNominalPos[0];
1551 initPos[1]=fNominalPos[1];
1552 initPos[2]=fNominalPos[2];
1553 }
1554
6d8df534 1555 Int_t nTrksSel = (Int_t)fTrkArraySel.GetEntries();
1556 if(nTrksSel==1) useConstraint=kTRUE;
3f2db92f 1557 AliDebug(1,Form("--- VertexFitter(): start (%d,%d,%d)",vfit,chiCalc,useWeights));
919e537f 1558 AliDebug(1,Form(" Number of tracks in array: %d\n",nTrksSel));
1559 AliDebug(1,Form(" Minimum # tracks required in fit: %d\n",fMinTracks));
1560 AliDebug(1,Form(" Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]));
1561 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 1562
0bf2e2b4 1563 // special treatment for few-tracks fits (e.g. secondary vertices)
3f2db92f 1564 Bool_t uUi3by3 = kFALSE; if(nTrksSel<5 && !useConstraint && !useWeights) uUi3by3 = kTRUE;
0bf2e2b4 1565
ec1be5d5 1566 Int_t i,j,k,step=0;
3f2db92f 1567 static TMatrixD rv(3,1);
1568 static TMatrixD vV(3,3);
ec1be5d5 1569 rv(0,0) = initPos[0];
1570 rv(1,0) = initPos[1];
3f2db92f 1571 rv(2,0) = initPos[2];
ec1be5d5 1572 Double_t xlStart,alpha;
3f2db92f 1573 Int_t nTrksUsed = 0;
1574 Double_t chi2=0,chi2i,chi2b;
6d8df534 1575 AliExternalTrackParam *t = 0;
ec1be5d5 1576 Int_t failed = 0;
1577
50ff0bcd 1578 // initial vertex covariance matrix
1579 TMatrixD vVb(3,3);
1580 vVb(0,0) = fNominalCov[0];
1581 vVb(0,1) = fNominalCov[1];
1582 vVb(0,2) = 0.;
1583 vVb(1,0) = fNominalCov[1];
1584 vVb(1,1) = fNominalCov[2];
1585 vVb(1,2) = 0.;
1586 vVb(2,0) = 0.;
1587 vVb(2,1) = 0.;
1588 vVb(2,2) = fNominalCov[5];
1589 TMatrixD vVbInv(TMatrixD::kInverted,vVb);
1590 TMatrixD rb(3,1);
1591 rb(0,0) = fNominalPos[0];
1592 rb(1,0) = fNominalPos[1];
1593 rb(2,0) = fNominalPos[2];
1594 TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb);
3f2db92f 1595 //
1596 int currBC = fVert.GetBC();
1597 //
eae1677e 1598 // 2 steps:
1599 // 1st - estimate of vtx using all tracks
1600 // 2nd - estimate of global chi2
3f2db92f 1601 //
eae1677e 1602 for(step=0; step<2; step++) {
3f2db92f 1603 if (step==0 && !vfit) continue;
1604 else if (step==1 && !chiCalc) continue;
ec1be5d5 1605 chi2 = 0.;
6d8df534 1606 nTrksUsed = 0;
3f2db92f 1607 fMVWSum = fMVWE2 = 0;
1608 if(step==1) { initPos[0]=rv(0,0); initPos[1]=rv(1,0); initPos[2]=rv(2,0);}
1609 AliDebug(2,Form("Step%d: inipos: %+f %+f %+f MinTr: %d, Sig2:%.2f)",step,initPos[0],initPos[1],initPos[2],fMinTracks,fMVSigma2));
1610 //
ec1be5d5 1611 TMatrixD sumWiri(3,1);
1612 TMatrixD sumWi(3,3);
1613 for(i=0; i<3; i++) {
1614 sumWiri(i,0) = 0.;
1615 for(j=0; j<3; j++) sumWi(j,i) = 0.;
1616 }
1617
fee9ef0d 1618 // mean vertex constraint
1619 if(useConstraint) {
50ff0bcd 1620 for(i=0;i<3;i++) {
1621 sumWiri(i,0) += vVbInvrb(i,0);
1622 for(k=0;k<3;k++) sumWi(i,k) += vVbInv(i,k);
07680cae 1623 }
fee9ef0d 1624 // chi2
1625 TMatrixD deltar = rv; deltar -= rb;
1626 TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar);
1627 chi2b = deltar(0,0)*vVbInvdeltar(0,0)+
1628 deltar(1,0)*vVbInvdeltar(1,0)+
1629 deltar(2,0)*vVbInvdeltar(2,0);
1630 chi2 += chi2b;
07680cae 1631 }
1632
ec1be5d5 1633 // loop on tracks
6d8df534 1634 for(k=0; k<nTrksSel; k++) {
3f2db92f 1635 //
ec1be5d5 1636 // get track from track array
6d8df534 1637 t = (AliExternalTrackParam*)fTrkArraySel.At(k);
3f2db92f 1638 if (useWeights && t->TestBit(kBitUsed)) continue;
1639 //
1640 int tBC = int(t->GetUniqueID()) - kTOFBCShift; // BC assigned to this track
1641 if (fSelectOnTOFBunchCrossing) {
1642 if (!fKeepAlsoUnflaggedTOFBunchCrossing) continue; // don't consider tracks with undefined BC
1643 if (currBC!=AliVTrack::kTOFBCNA && tBC!=AliVTrack::kTOFBCNA && tBC!=currBC) continue; // track does not match to current BCid
1644 }
ec1be5d5 1645 alpha = t->GetAlpha();
1646 xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
fee9ef0d 1647 // to vtxSeed (from finder)
6d8df534 1648 t->PropagateTo(xlStart,GetFieldkG());
fee9ef0d 1649
ec1be5d5 1650 // vector of track global coordinates
1651 TMatrixD ri(3,1);
fee9ef0d 1652 // covariance matrix of ri
1653 TMatrixD wWi(3,3);
1654
1655 // get space point from track
0bf2e2b4 1656 if(!TrackToPoint(t,ri,wWi,uUi3by3)) continue;
ec1be5d5 1657
1658 // track chi2
1659 TMatrixD deltar = rv; deltar -= ri;
1660 TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
1661 chi2i = deltar(0,0)*wWideltar(0,0)+
1662 deltar(1,0)*wWideltar(1,0)+
1663 deltar(2,0)*wWideltar(2,0);
3f2db92f 1664 //
1665 if (useWeights) {
1666 //double sg = TMath::Sqrt(fMVSigma2);
1667 //double chi2iw = (deltar(0,0)*wWideltar(0,0)+deltar(1,0)*wWideltar(1,0))/sg + deltar(2,0)*wWideltar(2,0)/fMVSigma2;
1668 //double wgh = (1-chi2iw/fMVTukey2);
1669 double chi2iw = chi2i;
1670 double wgh = (1-chi2iw/fMVTukey2/fMVSigma2);
1671
1672 if (wgh<kTiny) wgh = 0;
1673 else if (useWeights==2) wgh = 1.; // use as binary weight
1674 if (step==1) ((AliESDtrack*)t)->SetBit(kBitUsed, wgh>0);
1675 if (wgh<kTiny) continue; // discard the track
1676 wWi *= wgh; // RS: use weight?
1677 fMVWSum += wgh;
1678 fMVWE2 += wgh*chi2iw;
1679 }
ec1be5d5 1680 // add to total chi2
3f2db92f 1681 if (fSelectOnTOFBunchCrossing && tBC!=AliVTrack::kTOFBCNA && currBC<0) currBC = tBC;
1682 //
ec1be5d5 1683 chi2 += chi2i;
3f2db92f 1684 TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
ec1be5d5 1685 sumWiri += wWiri;
1686 sumWi += wWi;
1687
6d8df534 1688 nTrksUsed++;
ec1be5d5 1689 } // end loop on tracks
1690
6d8df534 1691 if(nTrksUsed < fMinTracks) {
ec1be5d5 1692 failed=1;
1693 continue;
1694 }
1695
1696 Double_t determinant = sumWi.Determinant();
f09c879d 1697 if(determinant < fMinDetFitter) {
919e537f 1698 AliDebug(1,Form("det(V) = %f (<%f)\n",determinant,fMinDetFitter));
ec1be5d5 1699 failed=1;
1700 continue;
1701 }
1702
0bf2e2b4 1703 if(step==0) {
1704 // inverted of weights matrix
1705 TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1706 vV = invsumWi;
1707 // position of primary vertex
1708 rv.Mult(vV,sumWiri);
1709 }
eae1677e 1710 } // end loop on the 2 steps
ec1be5d5 1711
ec1be5d5 1712 if(failed) {
3f2db92f 1713 fVert.SetNContributors(-1);
1714 if (chiCalc) {
1715 TooFewTracks();
1716 if (fCurrentVertex) fVert = *fCurrentVertex; // RS
1717 }
ec1be5d5 1718 return;
1719 }
3f2db92f 1720 //
ec1be5d5 1721 Double_t position[3];
1722 position[0] = rv(0,0);
1723 position[1] = rv(1,0);
1724 position[2] = rv(2,0);
1725 Double_t covmatrix[6];
1726 covmatrix[0] = vV(0,0);
1727 covmatrix[1] = vV(0,1);
1728 covmatrix[2] = vV(1,1);
1729 covmatrix[3] = vV(0,2);
1730 covmatrix[4] = vV(1,2);
1731 covmatrix[5] = vV(2,2);
1732
dc3719f3 1733 // for correct chi2/ndf, count constraint as additional "track"
6d8df534 1734 if(fConstraint) nTrksUsed++;
3f2db92f 1735 //
1736 if (vfit && !chiCalc) { // RS: special mode for multi-vertex finder
1737 fVert.SetXYZ(position);
1738 fVert.SetCovarianceMatrix(covmatrix);
1739 fVert.SetNContributors(nTrksUsed);
1740 return;
1741 }
ec1be5d5 1742 // store data in the vertex object
5b78b791 1743 if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
6d8df534 1744 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
3f2db92f 1745 fCurrentVertex->SetBC(currBC);
1746 fVert = *fCurrentVertex; // RS
919e537f 1747 AliDebug(1," Vertex after fit:");
1748 AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
1749 AliDebug(1,"--- VertexFitter(): finish\n");
1750
ec1be5d5 1751
1752 return;
1753}
50ff0bcd 1754//----------------------------------------------------------------------------
f5740e9a 1755AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(const TObjArray *trkArray,
6d8df534 1756 UShort_t *id,
1757 Bool_t optUseFitter,
97669588 1758 Bool_t optPropagate,
1759 Bool_t optUseDiamondConstraint)
fee9ef0d 1760{
eae1677e 1761//
6d8df534 1762// Return vertex from tracks (AliExternalTrackParam) in array
eae1677e 1763//
5b78b791 1764 fCurrentVertex = 0;
97669588 1765 // set optUseDiamondConstraint=TRUE only if you are reconstructing the
1766 // primary vertex!
1767 if(optUseDiamondConstraint) {
1768 SetConstraintOn();
1769 } else {
1770 SetConstraintOff();
1771 }
fee9ef0d 1772
50ff0bcd 1773 // get tracks and propagate them to initial vertex position
6d8df534 1774 fIdSel = new UShort_t[(Int_t)trkArray->GetEntriesFast()];
1775 Int_t nTrksSel = PrepareTracks(*trkArray,id,0);
97669588 1776 if((!optUseDiamondConstraint && nTrksSel<TMath::Max(2,fMinTracks)) ||
1777 (optUseDiamondConstraint && nTrksSel<1)) {
6d8df534 1778 TooFewTracks();
fee9ef0d 1779 return fCurrentVertex;
50ff0bcd 1780 }
1781
6d8df534 1782 // vertex finder
97669588 1783 if(nTrksSel==1) {
1784 AliDebug(1,"Just one track");
1785 OneTrackVertFinder();
1786 } else {
1787 switch (fAlgo) {
1788 case 1: StrLinVertexFinderMinDist(1); break;
1789 case 2: StrLinVertexFinderMinDist(0); break;
1790 case 3: HelixVertexFinder(); break;
1791 case 4: VertexFinder(1); break;
1792 case 5: VertexFinder(0); break;
1793 default: printf("Wrong algorithm\n"); break;
1794 }
fee9ef0d 1795 }
919e537f 1796 AliDebug(1," Vertex finding completed\n");
fee9ef0d 1797
1798 // vertex fitter
6d8df534 1799 if(optUseFitter) {
1800 VertexFitter();
1801 } else {
fee9ef0d 1802 Double_t position[3]={fVert.GetXv(),fVert.GetYv(),fVert.GetZv()};
146c29df 1803 Double_t covmatrix[6];
1804 fVert.GetCovMatrix(covmatrix);
fee9ef0d 1805 Double_t chi2=99999.;
6d8df534 1806 Int_t nTrksUsed=fVert.GetNContributors();
1807 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
fee9ef0d 1808 }
1809 fCurrentVertex->SetDispersion(fVert.GetDispersion());
1810
1811
1812 // set indices of used tracks and propagate track to found vertex
50ff0bcd 1813 UShort_t *indices = 0;
6d8df534 1814 Double_t d0z0[2],covd0z0[3];
1815 AliExternalTrackParam *t = 0;
fee9ef0d 1816 if(fCurrentVertex->GetNContributors()>0) {
44a76311 1817 indices = new UShort_t[fTrkArraySel.GetEntriesFast()];
6d8df534 1818 for(Int_t jj=0; jj<(Int_t)fTrkArraySel.GetEntriesFast(); jj++) {
1819 indices[jj] = fIdSel[jj];
1820 t = (AliExternalTrackParam*)fTrkArraySel.At(jj);
1821 if(optPropagate && optUseFitter) {
fee9ef0d 1822 if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
6d8df534 1823 t->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
919e537f 1824 AliDebug(1,Form("Track %d propagated to found vertex",jj));
6d8df534 1825 } else {
fee9ef0d 1826 AliWarning("Found vertex outside beam pipe!");
1827 }
1828 }
50ff0bcd 1829 }
fee9ef0d 1830 fCurrentVertex->SetIndices(fCurrentVertex->GetNContributors(),indices);
50ff0bcd 1831 }
6d8df534 1832
1833 // clean up
4ce766eb 1834 if (indices) {delete [] indices; indices=NULL;}
6fefa5ce 1835 delete [] fIdSel; fIdSel=NULL;
6d8df534 1836 fTrkArraySel.Delete();
50ff0bcd 1837
fee9ef0d 1838 return fCurrentVertex;
eae1677e 1839}
3f2db92f 1840
50ff0bcd 1841//----------------------------------------------------------------------------
6d8df534 1842AliESDVertex* AliVertexerTracks::VertexForSelectedESDTracks(TObjArray *trkArray,
3f2db92f 1843 Bool_t optUseFitter,
97669588 1844 Bool_t optPropagate,
1845 Bool_t optUseDiamondConstraint)
1846
fee9ef0d 1847{
50ff0bcd 1848//
6d8df534 1849// Return vertex from array of ESD tracks
50ff0bcd 1850//
5834e9a1 1851
6d8df534 1852 Int_t nTrks = (Int_t)trkArray->GetEntriesFast();
1853 UShort_t *id = new UShort_t[nTrks];
1854
1855 AliESDtrack *esdt = 0;
50ff0bcd 1856 for(Int_t i=0; i<nTrks; i++){
6d8df534 1857 esdt = (AliESDtrack*)trkArray->At(i);
1858 id[i] = (UShort_t)esdt->GetID();
50ff0bcd 1859 }
1860
97669588 1861 VertexForSelectedTracks(trkArray,id,optUseFitter,optPropagate,optUseDiamondConstraint);
6d8df534 1862
4ce766eb 1863 delete [] id; id=NULL;
6d8df534 1864
1865 return fCurrentVertex;
50ff0bcd 1866}
3f2db92f 1867
1868//______________________________________________________
1869Bool_t AliVertexerTracks::FindNextVertexMV()
1870{
1871 // try to find a new vertex
1872 fMVSigma2 = fMVSig2Ini;
1873 double prevSig2 = fMVSigma2;
1874 double minDst2 = fMVMinDst*fMVMinDst;
1875 const double kSigLimit = 1.0;
1876 const double kSigLimitE = kSigLimit+1e-6;
1877 const double kPushFactor = 0.5;
1878 const int kMaxIter = 20;
1879 double push = kPushFactor;
1880 //
1881 int iter = 0;
1882 double posP[3]={0,0,0},pos[3]={0,0,0};
1883 fVert.GetXYZ(posP);
1884 //
1885
1886 do {
1887 fVert.SetBC(AliVTrack::kTOFBCNA);
1888 VertexFitter(kTRUE,kFALSE,1);
1889 if (fVert.GetNContributors()<fMinTracks) {
1890 AliDebug(3,Form("Failed in iteration %d: No Contributirs",iter));
1891 break;
1892 } // failed
1893 if (fMVWSum>0) fMVSigma2 = TMath::Max(fMVWE2/fMVWSum,kSigLimit);
1894 else {
1895 AliDebug(3,Form("Failed %d, no weithgs",iter));
1896 iter = kMaxIter+1;
1897 break;
1898 } // failed
1899 //
1900 double sigRed = (prevSig2-fMVSigma2)/prevSig2; // reduction of sigma2
1901 //
1902 fVert.GetXYZ(pos);
1903 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]);
1904 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));
1905 if ( (++iter<kMaxIter) && (sigRed<0 || sigRed<fMVMinSig2Red) && fMVSigma2>fMVMaxSigma2) {
1906 fMVSigma2 *= push; // stuck, push little bit
1907 push *= kPushFactor;
1908 if (fMVSigma2<1.) fMVSigma2 = 1.;
1909 AliDebug(3,Form("Pushed sigma2 to %f",fMVSigma2));
1910 }
1911 else if (dst2<minDst2 && ((sigRed<0 || sigRed<fMVMinSig2Red) || fMVSigma2<kSigLimitE)) break;
1912 //
1913 fVert.GetXYZ(posP); // fetch previous vertex position
1914 prevSig2 = fMVSigma2;
1915 } while(iter<kMaxIter);
1916 //
1917 if (fVert.GetNContributors()<0 || iter>kMaxIter || fMVSigma2>fMVMaxSigma2) {
1918 return kFALSE;
1919 }
1920 else {
1921 VertexFitter(kFALSE,kTRUE,fMVFinalWBinary ? 2:1); // final chi2 calculation
1922 int nv = fMVVertices->GetEntries();
1923 // create indices
1924 int ntrk = fTrkArraySel.GetEntries();
1925 int nindices = fCurrentVertex->GetNContributors() - (fConstraint ? 1:0);
1926 UShort_t *indices = 0;
1927 if (nindices>0) indices = new UShort_t[nindices];
1928 int nadded = 0;
1929 for (int itr=0;itr<ntrk;itr++) {
1930 AliExternalTrackParam* t = (AliExternalTrackParam*)fTrkArraySel[itr];
1931 if (t->TestBit(kBitAccounted) || !t->TestBit(kBitUsed)) continue; // already belongs to some vertex
1932 t->SetBit(kBitAccounted);
1933 indices[nadded++] = fIdSel[itr];
1934 }
1935 if (nadded!=nindices) printf("Mismatch : NInd: %d Nadd: %d\n",nindices,nadded);
1936 fCurrentVertex->SetIndices(nindices,indices);
1937 // set vertex title
bac16cc1 1938 TString title="VertexerTracksMVNoConstraint";
1939 if(fConstraint) title="VertexerTracksMVWithConstraint";
3f2db92f 1940 fCurrentVertex->SetTitle(title.Data());
1941 fMVVertices->AddLast(fCurrentVertex);
1942 AliDebug(3,Form("Added new vertex #%d NCont:%d XYZ: %f %f %f",nindices,nv,fCurrentVertex->GetX(),fCurrentVertex->GetY(),fCurrentVertex->GetZ()));
1943 if (indices) delete[] indices;
1944 fCurrentVertex = 0; // already attached to fMVVertices
1945 return kTRUE;
1946 }
1947}
1948
1949//______________________________________________________
1950void AliVertexerTracks::FindVerticesMV()
1951{
1952 // find and fit multiple vertices
1953 //
1954 double step = fMVScanStep>1 ? fMVScanStep : 1.;
1955 double zmx = 3*TMath::Sqrt(fNominalCov[5]);
1956 double zmn = -zmx;
12f9872e 1957 int nz = TMath::Nint((zmx-zmn)/step); if (nz<1) nz=1;
3f2db92f 1958 double dz = (zmx-zmn)/nz;
1959 int izStart=0;
1960 AliDebug(2,Form("%d seeds between %f and %f",nz,zmn+dz/2,zmx+dz/2));
1961 //
1962 if (!fMVVertices) fMVVertices = new TObjArray(10);
1963 fMVVertices->Clear();
1964 //
1965 int ntrLeft = (Int_t)fTrkArraySel.GetEntries();
1966 //
1967 double sig2Scan = fMVSig2Ini;
1968 Bool_t runMore = kTRUE;
1969 int cntWide = 0;
1970 while (runMore) {
1971 fMVSig2Ini = sig2Scan*1e3; // try wide search
1972 Bool_t found = kFALSE;
1973 cntWide++;
1974 fVert.SetNContributors(-1);
1975 fVert.SetXYZ(fNominalPos);
1976 AliDebug(3,Form("Wide search #%d Z= %f Sigma2=%f",cntWide,fNominalPos[2],fMVSig2Ini));
1977 if (FindNextVertexMV()) { // are there tracks left to consider?
1978 AliESDVertex* vtLast = (AliESDVertex*)fMVVertices->Last();
1979 if (vtLast && vtLast->GetNContributors()>0) ntrLeft -= vtLast->GetNContributors()-(fConstraint ? 1:0);
1980 if (ntrLeft<1) runMore = kFALSE;
1981 found = kTRUE;
1982 continue;
1983 }
1984 // if nothing is found, do narrow sig2ini scan
1985 fMVSig2Ini = sig2Scan;
1986 for (;izStart<nz;izStart++) {
1987 double zSeed = zmn+dz*(izStart+0.5);
1988 AliDebug(3,Form("Seed %d: Z= %f Sigma2=%f",izStart,zSeed,fMVSig2Ini));
1989 fVert.SetNContributors(-1);
1990 fVert.SetXYZ(fNominalPos);
1991 fVert.SetZv(zSeed);
1992 if (FindNextVertexMV()) { // are there tracks left to consider?
1993 AliESDVertex* vtLast = (AliESDVertex*)fMVVertices->Last();
1994 if (vtLast && vtLast->GetNContributors()>0) ntrLeft -= vtLast->GetNContributors()-(fConstraint ? 1:0);
1995 if (ntrLeft<1) runMore = kFALSE;
1996 found = kTRUE;
1997 break;
1998 }
1999 }
2000 runMore = found; // if nothing was found, no need for new iteration
2001 }
2002 fMVSig2Ini = sig2Scan;
2003 int nvFound = fMVVertices->GetEntriesFast();
2004 AliDebug(2,Form("Number of found vertices: %d",nvFound));
2005 if (nvFound<1) TooFewTracks();
2006 if (AliLog::GetGlobalDebugLevel()>0) fMVVertices->Print();
2007 //
2008}
2009
2010//______________________________________________________
2011void AliVertexerTracks::AnalyzePileUp(AliESDEvent* esdEv)
2012{
2013 // if multiple vertices are found, try to find the primary one and attach it as fCurrentVertex
2014 // then attach pile-up vertices directly to esdEv
2015 //
2016 int nFND = (fMVVertices && fMVVertices->GetEntriesFast()) ? fMVVertices->GetEntriesFast() : 0;
2017 if (nFND<1) { if (!fCurrentVertex) TooFewTracks(); return;} // no multiple vertices
2018 //
2019 int indCont[nFND];
2020 int nIndx[nFND];
2021 for (int iv=0;iv<nFND;iv++) {
2022 AliESDVertex* fnd = (AliESDVertex*)fMVVertices->At(iv);
2023 indCont[iv] = iv;
2024 nIndx[iv] = fnd->GetNIndices();
2025 }
2026 TMath::Sort(nFND, nIndx, indCont, kTRUE); // sort in decreasing order of Nindices
2027 double dists[nFND];
2028 int distOrd[nFND],indx[nFND];
2029 for (int iv=0;iv<nFND;iv++) {
2030 AliESDVertex* fndI = (AliESDVertex*)fMVVertices->At(indCont[iv]);
2031 if (fndI->GetStatus()<1) continue; // discarded
2032 int ncomp = 0;
2033 for (int jv=iv+1;jv<nFND;jv++) {
2034 AliESDVertex* fndJ = (AliESDVertex*)fMVVertices->At(indCont[jv]);
2035 if (fndJ->GetStatus()<1) continue;
2036 dists[ncomp] = fndI->GetWDist(fndJ)*fndJ->GetNIndices();
2037 distOrd[ncomp] = indCont[jv];
2038 indx[ncomp] = ncomp;
2039 ncomp++;
2040 }
2041 if (ncomp<1) continue;
2042 TMath::Sort(ncomp, dists, indx, kFALSE); // sort in increasing distance order
2043 for (int jv0=0;jv0<ncomp;jv0++) {
2044 int jv = distOrd[indx[jv0]];
2045 AliESDVertex* fndJ = (AliESDVertex*)fMVVertices->At(jv);
2046 if (dists[indx[jv0]]<fMVMaxWghNtr) { // candidate for split vertex
2047 //before eliminating the small close vertex, check if we should transfere its BC to largest one
2048 if (fndJ->GetBC()!=AliVTrack::kTOFBCNA && fndI->GetBC()==AliVTrack::kTOFBCNA) fndI->SetBC(fndJ->GetBC());
2049 //
2050 // leave the vertex only if both vertices have definit but different BC's, then this is not splitting.
2051 if ( (fndJ->GetBC()==fndI->GetBC()) || (fndJ->GetBC()==AliVTrack::kTOFBCNA)) fndJ->SetNContributors(-fndJ->GetNContributors());
2052 }
2053 }
2054 }
2055 //
2056 // select as a primary the largest vertex with BC=0, or the largest with BC non-ID
2057 int primBC0=-1,primNoBC=-1;
2058 for (int iv0=0;iv0<nFND;iv0++) {
2059 int iv = indCont[iv0];
2060 AliESDVertex* fndI = (AliESDVertex*)fMVVertices->At(iv);
2061 if (!fndI) continue;
2062 if (fndI->GetStatus()<1) {fMVVertices->RemoveAt(iv); delete fndI; continue;}
2063 if (primBC0<0 && fndI->GetBC()==0) primBC0 = iv;
2064 if (primNoBC<0 && fndI->GetBC()==AliVTrack::kTOFBCNA) primNoBC = iv;
2065 }
2066 //
2067 if (primBC0>=0) fCurrentVertex = (AliESDVertex*)fMVVertices->At(primBC0);
2068 if (!fCurrentVertex && primNoBC>=0) fCurrentVertex = (AliESDVertex*)fMVVertices->At(primNoBC);
2069 if (fCurrentVertex) fMVVertices->Remove(fCurrentVertex);
2070 else { // all vertices have BC>0, no primary vertex
2071 fCurrentVertex = new AliESDVertex();
2072 fCurrentVertex->SetNContributors(-3);
2073 fCurrentVertex->SetBC(AliVTrack::kTOFBCNA);
2074 }
2075 fCurrentVertex->SetID(-1);
2076 //
2077 // add pileup vertices
2078 int nadd = 0;
2079 for (int iv0=0;iv0<nFND;iv0++) {
2080 int iv = indCont[iv0];
2081 AliESDVertex* fndI = (AliESDVertex*)fMVVertices->At(iv);
2082 if (!fndI) continue;
2083 fndI->SetID(++nadd);
2084 esdEv->AddPileupVertexTracks(fndI);
2085 }
2086 //
2087 fMVVertices->Delete();
2088 //
2089}
2090
7878b3ae 2091//______________________________________________________
2092void AliVertexerTracks::FindAllVertices(Int_t nTrksOrig,
2093 const TObjArray *trkArrayOrig,
2094 Double_t* zTr,
2095 Double_t* err2zTr,
2096 UShort_t* idOrig){
2097
2098 // clusterize tracks using z coordinates of intersection with beam axis
2099 // and compute all vertices
2100
2101 UShort_t* posOrig=new UShort_t[nTrksOrig];
2102 for(Int_t iTr=0; iTr<nTrksOrig; iTr++) posOrig[iTr]=iTr;
2103
2104
2105 // sort points along Z
2106 AliDebug(1,Form("Sort points along Z, used tracks %d",nTrksOrig));
2107 for(Int_t iTr1=0; iTr1<nTrksOrig; iTr1++){
2108 for(Int_t iTr2=iTr1+1; iTr2<nTrksOrig; iTr2++){
2109 if(zTr[iTr1]>zTr[iTr2]){
2110 Double_t tmpz=zTr[iTr2];
2111 Double_t tmperr=err2zTr[iTr2];
2112 UShort_t tmppos=posOrig[iTr2];
2113 UShort_t tmpid=idOrig[iTr2];
2114 zTr[iTr2]=zTr[iTr1];
2115 err2zTr[iTr2]=err2zTr[iTr1];
2116 posOrig[iTr2]=posOrig[iTr1];
2117 idOrig[iTr2]=idOrig[iTr1];
2118 zTr[iTr1]=tmpz;
2119 err2zTr[iTr1]=tmperr;
2120 idOrig[iTr1]=tmpid;
2121 posOrig[iTr1]=tmppos;
2122 }
2123 }
2124 }
2125
2126 // clusterize
2127 Int_t nClusters=0;
2128 Int_t* firstTr=new Int_t[nTrksOrig];
2129 Int_t* lastTr=new Int_t[nTrksOrig];
2130
2131 firstTr[0]=0;
2132 for(Int_t iTr=0; iTr<nTrksOrig-1; iTr++){
2133 Double_t distz=zTr[iTr+1]-zTr[iTr];
2134 Double_t errdistz=TMath::Sqrt(err2zTr[iTr+1]+err2zTr[iTr]);
2135 if(errdistz<=0.000001) errdistz=0.000001;
2136 if(distz>fDeltaZCutForCluster || (distz/errdistz)>fnSigmaZCutForCluster){
2137 lastTr[nClusters]=iTr;
2138 firstTr[nClusters+1]=iTr+1;
2139 nClusters++;
2140 }
2141 }
2142 lastTr[nClusters]=nTrksOrig-1;
2143
2144 // Compute vertices
2145 AliDebug(1,Form("Number of found clusters %d",nClusters+1));
2146 Int_t nFoundVertices=0;
2147
2148 if (!fMVVertices) fMVVertices = new TObjArray(nClusters+1);
2149
2150 fMVVertices->Clear();
2151 TObjArray cluTrackArr(nTrksOrig);
2152 UShort_t *idTrkClu=new UShort_t[nTrksOrig];
2153 // Int_t maxContr=0;
2154 // Int_t maxPos=-1;
2155
2156 for(Int_t iClu=0; iClu<=nClusters; iClu++){
2157 Int_t nCluTracks=lastTr[iClu]-firstTr[iClu]+1;
2158 cluTrackArr.Clear();
2159 AliDebug(1,Form(" Vertex #%d tracks %d first tr %d last track %d",iClu,nCluTracks,firstTr[iClu],lastTr[iClu]));
2160 Int_t nSelTr=0;
2161 for(Int_t iTr=firstTr[iClu]; iTr<=lastTr[iClu]; iTr++){
2162 AliExternalTrackParam* t=(AliExternalTrackParam*)trkArrayOrig->At(posOrig[iTr]);
2163 if(TMath::Abs(t->GetZ()-zTr[iTr])>0.00001){
2164 AliError(Form("Clu %d Track %d zTrack=%f zVec=%f\n",iClu,iTr,t->GetZ(),zTr[iTr]));
2165 }
2166 cluTrackArr.AddAt(t,nSelTr);
2167 idTrkClu[nSelTr]=idOrig[iTr];
2168 AliDebug(1,Form(" Add track %d: id %d, z=%f",iTr,idOrig[iTr],zTr[iTr]));
2169 nSelTr++;
2170 }
2171 AliESDVertex* vert=FindPrimaryVertex(&cluTrackArr,idTrkClu);
2172 AliDebug(1,Form("Found vertex in z=%f with %d contributors",vert->GetZv(),
2173 vert->GetNContributors()));
2174
2175 fCurrentVertex=0;
2176 if(vert->GetNContributors()>0){
2177 nFoundVertices++;
2178 fMVVertices->AddLast(vert);
2179 }
2180 // if(vert->GetNContributors()>maxContr){
2181 // maxContr=vert->GetNContributors();
2182 // maxPos=nFoundVertices-1;
2183 // }
2184 }
2185
2186 AliDebug(1,Form("Number of found vertices %d (%d)",nFoundVertices,fMVVertices->GetEntriesFast()));
2187 // if(maxPos>=0 && maxContr>0){
2188 // AliESDVertex* vtxMax=(AliESDVertex*)fMVVertices->At(maxPos);
2189 // if(fCurrentVertex) delete fCurrentVertex;
2190 // fCurrentVertex=new AliESDVertex(*vtxMax);
2191 // }
2192
2193 delete [] firstTr;
2194 delete [] lastTr;
2195 delete [] idTrkClu;
2196 delete [] posOrig;
2197
2198 return;
2199
2200}