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