Add tpc seed to the list of object for calibration in the esd track (M.Ivanov)
[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//-----------------------------------------------------------------
17// Implementation of the vertexer from ESD tracks
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:
25// F.Prino, Torino, prino@to.infn.it
26//-----------------------------------------------------------------
27
28//---- Root headers --------
29#include <TTree.h>
30//---- AliRoot headers -----
31#include "AliStrLine.h"
32#include "AliVertexerTracks.h"
33#include "AliESDtrack.h"
34
35ClassImp(AliVertexerTracks)
36
37
38//----------------------------------------------------------------------------
39AliVertexerTracks::AliVertexerTracks():
40 TObject(),fVert()
41{
42//
43// Default constructor
44//
45 SetVtxStart();
46 SetMinTracks();
47 fDCAcut=0;
48 fAlgo=1;
49}
50//-----------------------------------------------------------------------------
51AliVertexerTracks::AliVertexerTracks(Double_t xStart, Double_t yStart):
52 TObject(),fVert()
53{
54//
55// Standard constructor
56//
57 SetVtxStart(xStart,yStart);
58 SetMinTracks();
59 fDCAcut=0;
60 fAlgo=1;
61}
62//-----------------------------------------------------------------------------
63AliVertexerTracks::~AliVertexerTracks() {
64 // Default Destructor
65 // The objects poited by the following pointers are not owned
66 // by this class and are not deleted
67}
68//----------------------------------------------------------------------------
69Int_t AliVertexerTracks::PrepareTracks(TTree &trkTree) {
70//
71// Propagate tracks to initial vertex position and store them in a TObjArray
72//
73 Double_t alpha,xlStart;
74 Int_t nTrks = 0;
75
76 Int_t nEntries = (Int_t)trkTree.GetEntries();
77 if(!fTrkArray.IsEmpty()) fTrkArray.Clear();
78 fTrkArray.Expand(nEntries);
79
80 for(Int_t i=0; i<nEntries; i++) {
81 AliESDtrack *track = new AliESDtrack;
82 trkTree.SetBranchAddress("tracks",&track);
83 trkTree.GetEvent(i);
84
85 // propagate track to vtxSeed
86 alpha = track->GetAlpha();
87 xlStart = fNominalPos[0]*TMath::Cos(alpha)+fNominalPos[1]*TMath::Sin(alpha);
88 track->PropagateTo(xlStart,GetField()); // to vtxSeed
89 fTrkArray.AddLast(track);
90 nTrks++;
91 }
92 return nTrks;
93}
94//----------------------------------------------------------------------------
95AliVertex* AliVertexerTracks::VertexForSelectedTracks(TTree *trkTree) {
96//
97// Return vertex from tracks in trkTree
98//
99
100 // get tracks and propagate them to initial vertex position
101 Int_t nTrks = PrepareTracks(*trkTree);
102 if(nTrks < fMinTracks) {
103 printf("TooFewTracks\n");
104 Double_t vtx[3]={0,0,0};
105 fVert.SetXYZ(vtx);
106 fVert.SetDispersion(999);
107 fVert.SetNContributors(-5);
108 return &fVert;
109 }
110
111 // Set initial vertex position from ESD
112 if(fAlgo==1) StrLinVertexFinderMinDist(1);
113 if(fAlgo==2) StrLinVertexFinderMinDist(0);
114 if(fAlgo==3) HelixVertexFinder();
115 if(fAlgo==4) VertexFinder(1);
116 if(fAlgo==5) VertexFinder(0);
117 return &fVert;
118}
119//----------------------------------------------------------------------------
120AliVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray) {
121//
122// Return vertex from array of tracks
123//
124
125 // get tracks and propagate them to initial vertex position
126 Int_t nTrks = trkArray->GetEntriesFast();
127 if(nTrks < fMinTracks) {
128 printf("TooFewTracks\n");
129 Double_t vtx[3]={0,0,0};
130 fVert.SetXYZ(vtx);
131 fVert.SetDispersion(999);
132 fVert.SetNContributors(-5);
133 return &fVert;
134 }
135 TTree *trkTree = new TTree("TreeT","tracks");
136 AliESDtrack *esdTrack = 0;
137 trkTree->Branch("tracks","AliESDtrack",&esdTrack);
138 for(Int_t i=0; i<nTrks; i++){
139 esdTrack = (AliESDtrack*)trkArray->At(i);
140 trkTree->Fill();
141 }
142
2e5c8fc7 143 AliVertex *vtx = VertexForSelectedTracks(trkTree);
144 delete trkTree;
145 return vtx;
2d57349e 146}
147
148//---------------------------------------------------------------------------
149void AliVertexerTracks::VertexFinder(Int_t optUseWeights) {
150
151 // Get estimate of vertex position in (x,y) from tracks DCA
152
153 Double_t initPos[3];
154 initPos[2] = 0.;
155 for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
156 Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
157 Double_t aver[3]={0.,0.,0.};
158 Double_t aversq[3]={0.,0.,0.};
159 Double_t sigmasq[3]={0.,0.,0.};
160 Double_t sigma=0;
161 Int_t ncombi = 0;
162 AliESDtrack *track1;
163 AliESDtrack *track2;
164 Double_t alpha,mindist;
165 Double_t field=GetField();
166
167 for(Int_t i=0; i<nacc; i++){
168 track1 = (AliESDtrack*)fTrkArray.At(i);
169 alpha=track1->GetAlpha();
170 mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
171 AliStrLine *line1 = new AliStrLine();
172 track1->ApproximateHelixWithLine(mindist,field,line1);
173
174 for(Int_t j=i+1; j<nacc; j++){
175 track2 = (AliESDtrack*)fTrkArray.At(j);
176 alpha=track2->GetAlpha();
177 mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
178 AliStrLine *line2 = new AliStrLine();
179 track2->ApproximateHelixWithLine(mindist,field,line2);
180 Double_t distCA=line2->GetDCA(line1);
181 if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
182 Double_t pnt1[3],pnt2[3],crosspoint[3];
183
184 if(optUseWeights<=0){
185 Int_t retcode = line2->Cross(line1,crosspoint);
186 if(retcode>=0){
187 ncombi++;
188 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
189 for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
190 }
191 }
192 if(optUseWeights>0){
193 Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2);
194 if(retcode>=0){
195 Double_t alpha, cs, sn;
196 alpha=track1->GetAlpha();
197 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
198 Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
199 alpha=track2->GetAlpha();
200 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
201 Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
202 Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
203 Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
204 Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
205 Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
206 crosspoint[0]=wx1*pnt1[0] + wx2*pnt2[0];
207 crosspoint[1]=wy1*pnt1[1] + wy2*pnt2[1];
208 crosspoint[2]=wz1*pnt1[2] + wz2*pnt2[2];
209
210 ncombi++;
211 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
212 for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
213 }
214 }
215 }
216 delete line2;
217 }
218 delete line1;
219 }
220 if(ncombi>0){
221 for(Int_t jj=0;jj<3;jj++){
222 initPos[jj] = aver[jj]/ncombi;
223 aversq[jj]/=ncombi;
224 sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
225 sigma+=sigmasq[jj];
226 }
227 sigma=TMath::Sqrt(TMath::Abs(sigma));
228 }
229 else {
230 Warning("VertexFinder","Finder did not succed");
231 sigma=999;
232 }
233 fVert.SetXYZ(initPos);
234 fVert.SetDispersion(sigma);
235 fVert.SetNContributors(ncombi);
236}
237//---------------------------------------------------------------------------
238void AliVertexerTracks::HelixVertexFinder() {
239
240 // Get estimate of vertex position in (x,y) from tracks DCA
241
242
243 Double_t initPos[3];
244 initPos[2] = 0.;
245 for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
246 Double_t field=GetField();
247
248 Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
249
250 Double_t aver[3]={0.,0.,0.};
251 Double_t averquad[3]={0.,0.,0.};
252 Double_t sigmaquad[3]={0.,0.,0.};
253 Double_t sigma=0;
254 Int_t ncombi = 0;
255 AliESDtrack *track1;
256 AliESDtrack *track2;
257 Double_t distCA;
258 Double_t x, par[5];
259 Double_t alpha, cs, sn;
260 Double_t crosspoint[3];
261 for(Int_t i=0; i<nacc; i++){
262 track1 = (AliESDtrack*)fTrkArray.At(i);
263
264
265 for(Int_t j=i+1; j<nacc; j++){
266 track2 = (AliESDtrack*)fTrkArray.At(j);
267
268 distCA=track2->PropagateToDCA(track1,field);
269
270 if(fDCAcut<=0 ||(fDCAcut>0&&distCA<fDCAcut)){
271 track1->GetExternalParameters(x,par);
272 alpha=track1->GetAlpha();
273 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
274 Double_t x1=x*cs - par[0]*sn;
275 Double_t y1=x*sn + par[0]*cs;
276 Double_t z1=par[1];
277 Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
278 track2->GetExternalParameters(x,par);
279 alpha=track2->GetAlpha();
280 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
281 Double_t x2=x*cs - par[0]*sn;
282 Double_t y2=x*sn + par[0]*cs;
283 Double_t z2=par[1];
284 Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
285 Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
286 Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
287 Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
288 Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
289 crosspoint[0]=wx1*x1 + wx2*x2;
290 crosspoint[1]=wy1*y1 + wy2*y2;
291 crosspoint[2]=wz1*z1 + wz2*z2;
292
293 ncombi++;
294 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
295 for(Int_t jj=0;jj<3;jj++)averquad[jj]+=(crosspoint[jj]*crosspoint[jj]);
296 }
297 }
298
299 }
300 if(ncombi>0){
301 for(Int_t jj=0;jj<3;jj++){
302 initPos[jj] = aver[jj]/ncombi;
303 averquad[jj]/=ncombi;
304 sigmaquad[jj]=averquad[jj]-initPos[jj]*initPos[jj];
305 sigma+=sigmaquad[jj];
306 }
307 sigma=TMath::Sqrt(TMath::Abs(sigma));
308 }
309 else {
310 Warning("HelixVertexFinder","Finder did not succed");
311 sigma=999;
312 }
313 fVert.SetXYZ(initPos);
314 fVert.SetDispersion(sigma);
315 fVert.SetNContributors(ncombi);
316}
317//---------------------------------------------------------------------------
318void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights){
319
320 // Calculate the point at minimum distance to prepared tracks
321
322 Double_t initPos[3];
323 initPos[2] = 0.;
324 Double_t sigma=0;
325 for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
326 const Int_t knacc = (Int_t)fTrkArray.GetEntriesFast();
327 Double_t field=GetField();
328
329 AliESDtrack *track1;
330 Double_t (*vectP0)[3]=new Double_t [knacc][3];
331 Double_t (*vectP1)[3]=new Double_t [knacc][3];
332
333 Double_t sum[3][3];
334 Double_t dsum[3]={0,0,0};
335 for(Int_t i=0;i<3;i++)
336 for(Int_t j=0;j<3;j++)sum[i][j]=0;
337 for(Int_t i=0; i<knacc; i++){
338 track1 = (AliESDtrack*)fTrkArray.At(i);
339 Double_t alpha=track1->GetAlpha();
340 Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
341 AliStrLine *line1 = new AliStrLine();
342 track1->ApproximateHelixWithLine(mindist,field,line1);
343
344 Double_t p0[3],cd[3];
345 line1->GetP0(p0);
346 line1->GetCd(cd);
347 Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
348 vectP0[i][0]=p0[0];
349 vectP0[i][1]=p0[1];
350 vectP0[i][2]=p0[2];
351 vectP1[i][0]=p1[0];
352 vectP1[i][1]=p1[1];
353 vectP1[i][2]=p1[2];
354
355 Double_t matr[3][3];
356 Double_t dknow[3];
357 if(optUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
358 if(optUseWeights==1){
359 Double_t sigmasq[3];
360 sigmasq[0]=track1->GetSigmaY2();
361 sigmasq[1]=track1->GetSigmaY2();
362 sigmasq[2]=track1->GetSigmaZ2();
363 GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
364 }
365
366 for(Int_t iii=0;iii<3;iii++){
367 dsum[iii]+=dknow[iii];
368 for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
369 }
370 delete line1;
371 }
372
373 Double_t vett[3][3];
374 Double_t det=GetDeterminant3X3(sum);
375
376 if(det!=0){
377 for(Int_t zz=0;zz<3;zz++){
378 for(Int_t ww=0;ww<3;ww++){
379 for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
380 }
381 for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
382 initPos[zz]=GetDeterminant3X3(vett)/det;
383 }
384
385
386 for(Int_t i=0; i<knacc; i++){
387 Double_t p0[3]={0,0,0},p1[3]={0,0,0};
388 for(Int_t ii=0;ii<3;ii++){
389 p0[ii]=vectP0[i][ii];
390 p1[ii]=vectP1[i][ii];
391 }
392 sigma+=GetStrLinMinDist(p0,p1,initPos);
393 }
394
395 sigma=TMath::Sqrt(sigma);
396 }else{
397 Warning("StrLinVertexFinderMinDist","Finder did not succed");
398 sigma=999;
399 }
400 delete vectP0;
401 delete vectP1;
402 fVert.SetXYZ(initPos);
403 fVert.SetDispersion(sigma);
404 fVert.SetNContributors(knacc);
405}
406//_______________________________________________________________________
407Double_t AliVertexerTracks::GetDeterminant3X3(Double_t matr[][3]){
408 //
409 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];
410 return det;
411}
412//____________________________________________________________________________
413void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t m[][3],Double_t *d){
414
415 //
416 Double_t x12=p0[0]-p1[0];
417 Double_t y12=p0[1]-p1[1];
418 Double_t z12=p0[2]-p1[2];
419 Double_t kk=x12*x12+y12*y12+z12*z12;
420 m[0][0]=2-2/kk*x12*x12;
421 m[0][1]=-2/kk*x12*y12;
422 m[0][2]=-2/kk*x12*z12;
423 m[1][0]=-2/kk*x12*y12;
424 m[1][1]=2-2/kk*y12*y12;
425 m[1][2]=-2/kk*y12*z12;
426 m[2][0]=-2/kk*x12*z12;
427 m[2][1]=-2*y12*z12;
428 m[2][2]=2-2/kk*z12*z12;
429 d[0]=2*p0[0]-2/kk*p0[0]*x12*x12-2/kk*p0[2]*x12*z12-2/kk*p0[1]*x12*y12;
430 d[1]=2*p0[1]-2/kk*p0[1]*y12*y12-2/kk*p0[0]*x12*y12-2/kk*p0[2]*z12*y12;
431 d[2]=2*p0[2]-2/kk*p0[2]*z12*z12-2/kk*p0[0]*x12*z12-2/kk*p0[1]*z12*y12;
432
433}
434//____________________________________________________________________________
435void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t *sigmasq,Double_t m[][3],Double_t *d){
436 //
437 Double_t x12=p1[0]-p0[0];
438 Double_t y12=p1[1]-p0[1];
439 Double_t z12=p1[2]-p0[2];
440
441 Double_t den= x12*x12*sigmasq[1]*sigmasq[2]+y12*y12*sigmasq[0]*sigmasq[2]+z12*z12*sigmasq[0]*sigmasq[1];
442
443 Double_t kk= 2*(x12*x12/sigmasq[0]+y12*y12/sigmasq[1]+z12*z12/sigmasq[2]);
444
445 Double_t cc[3];
446 cc[0]=-x12/sigmasq[0];
447 cc[1]=-y12/sigmasq[1];
448 cc[2]=-z12/sigmasq[2];
449
450 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;
451
452 Double_t ss= -p0[0]*cc[0]-p0[1]*cc[1]-p0[2]*cc[2];
453
454 Double_t aa[3];
455 aa[0]=x12*sigmasq[1]*sigmasq[2]/den;
456 aa[1]=y12*sigmasq[0]*sigmasq[2]/den;
457 aa[2]=z12*sigmasq[0]*sigmasq[1]/den;
458
459 m[0][0]=aa[0]*(aa[0]*kk+2*cc[0])+2*cc[0]*aa[0]+2/sigmasq[0];
460 m[0][1]=aa[1]*(aa[0]*kk+2*cc[0])+2*cc[1]*aa[0];
461 m[0][2]=aa[2]*(aa[0]*kk+2*cc[0])+2*cc[2]*aa[0];
462
463 m[1][0]=aa[0]*(aa[1]*kk+2*cc[1])+2*cc[0]*aa[1];
464 m[1][1]=aa[1]*(aa[1]*kk+2*cc[1])+2*cc[1]*aa[1]+2/sigmasq[1];
465 m[1][2]=aa[2]*(aa[1]*kk+2*cc[1])+2*cc[2]*aa[1];
466
467 m[2][0]=aa[0]*(aa[2]*kk+2*cc[2])+2*cc[0]*aa[2];
468 m[2][1]=aa[1]*(aa[2]*kk+2*cc[2])+2*cc[1]*aa[2];
469 m[2][2]=aa[2]*(aa[2]*kk+2*cc[2])+2*cc[2]*aa[2]+2/sigmasq[2];
470
471 d[0]=-ww*(aa[0]*kk+2*cc[0])-2*ss*aa[0]+2*p0[0]/sigmasq[0];
472 d[1]=-ww*(aa[1]*kk+2*cc[1])-2*ss*aa[1]+2*p0[1]/sigmasq[1];
473 d[2]=-ww*(aa[2]*kk+2*cc[2])-2*ss*aa[2]+2*p0[2]/sigmasq[2];
474
475}
476//_____________________________________________________________________________
477Double_t AliVertexerTracks::GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t *x0){
478 //
479 Double_t x12=p0[0]-p1[0];
480 Double_t y12=p0[1]-p1[1];
481 Double_t z12=p0[2]-p1[2];
482 Double_t x10=p0[0]-x0[0];
483 Double_t y10=p0[1]-x0[1];
484 Double_t z10=p0[2]-x0[2];
485 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);
486}