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