Changes required by Effective C++
[u/mrichter/AliRoot.git] / STEER / AliVertexerTracks.cxx
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 #include <TMatrixD.h>
31 //---- AliRoot headers -----
32 #include "AliStrLine.h"
33 #include "AliVertexerTracks.h"
34 #include "AliESD.h"
35 #include "AliESDtrack.h"
36
37 ClassImp(AliVertexerTracks)
38
39
40 //----------------------------------------------------------------------------
41 AliVertexerTracks::AliVertexerTracks():
42 TObject(),
43 fVert(),
44 fCurrentVertex(0),
45 fMinTracks(2),
46 fMinITSClusters(5),
47 fTrkArray(),
48 fTrksToSkip(0),
49 fNTrksToSkip(0),
50 fDCAcut(0),
51 fAlgo(1),
52 fNSigma(3),
53 fDebug(0)
54 {
55 //
56 // Default constructor
57 //
58   SetVtxStart();
59   SetVtxStartSigma();
60   SetMinTracks();
61   SetMinITSClusters();
62   SetNSigmad0(); 
63 }
64 //-----------------------------------------------------------------------------
65 AliVertexerTracks::AliVertexerTracks(Double_t xStart, Double_t yStart):
66 TObject(),
67 fVert(),
68 fCurrentVertex(0),
69 fMinTracks(2),
70 fMinITSClusters(5),
71 fTrkArray(),
72 fTrksToSkip(0),
73 fNTrksToSkip(0),
74 fDCAcut(0),
75 fAlgo(1),
76 fNSigma(3),
77 fDebug(0)
78 {
79 //
80 // Standard constructor
81 //
82   SetVtxStart(xStart,yStart);
83   SetVtxStartSigma();
84   SetMinTracks();
85   SetMinITSClusters();
86   SetNSigmad0(); 
87 }
88 //-----------------------------------------------------------------------------
89 AliVertexerTracks::~AliVertexerTracks() {
90   // Default Destructor
91   // The objects pointed by the following pointer are not owned
92   // by this class and are not deleted
93   fCurrentVertex = 0;
94   delete[] fTrksToSkip;
95 }
96 //---------------------------------------------------------------------------
97 void AliVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) {
98 //
99 // Mark the tracks not ot be used in the vertex finding
100 //
101   fNTrksToSkip = n;
102   fTrksToSkip = new Int_t[n]; 
103   for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i]; 
104   return; 
105 }
106 //----------------------------------------------------------------------------
107 Int_t AliVertexerTracks::PrepareTracks(TTree &trkTree, Int_t OptImpParCut) {
108 //
109 // Propagate tracks to initial vertex position and store them in a TObjArray
110 //
111   Double_t maxd0rphi = 3.;  
112   Int_t    nTrks    = 0;
113   Double_t sigmaCurr[3],sigmaVtx[3],posVtx[3];
114   Double_t covVtx[6],covVtxXY;
115   Float_t  d0z0[2],covd0z0[3]; 
116   Double_t momentum[3],postrk[3];
117   Double_t pt,sigma,sigmavtxphi,phitrk;
118   Double_t field=GetField();
119
120   AliESDVertex *initVertex = new AliESDVertex(fNominalPos,fNominalSigma);
121
122   Int_t    nEntries = (Int_t)trkTree.GetEntries();
123   if(!fTrkArray.IsEmpty()) fTrkArray.Clear();
124   fTrkArray.Expand(nEntries);
125
126   if(fDebug) {
127     printf(" PrepareTracks()\n");
128   }
129
130   for(Int_t i=0; i<nEntries; i++) {
131     AliESDtrack *track = new AliESDtrack; 
132     trkTree.SetBranchAddress("tracks",&track);
133     trkTree.GetEvent(i);
134
135
136
137     // propagate track to vertex
138     if(OptImpParCut==1) { // OptImpParCut==1
139       track->RelateToVertex(initVertex,field,100.);
140       initVertex->GetXYZ(posVtx);
141       initVertex->GetSigmaXYZ(sigmaVtx);
142       covVtxXY = 0.;
143     } else {              // OptImpParCut==2
144       fCurrentVertex->GetSigmaXYZ(sigmaCurr);
145       if((sigmaCurr[0]+sigmaCurr[1])<(fNominalSigma[0]+fNominalSigma[1])) {
146         track->RelateToVertex(fCurrentVertex,field,100.);
147         fCurrentVertex->GetXYZ(posVtx);
148         fCurrentVertex->GetSigmaXYZ(sigmaVtx);
149         fCurrentVertex->GetCovMatrix(covVtx);
150         covVtxXY = covVtx[1];
151       } else {
152         track->RelateToVertex(initVertex,field,100.);
153         initVertex->GetXYZ(posVtx);
154         initVertex->GetSigmaXYZ(sigmaVtx);
155         covVtxXY = 0.;
156       }
157     }
158
159     // select tracks with d0rphi < maxd0rphi
160     track->GetImpactParameters(d0z0,covd0z0);
161
162     track->GetXYZ(postrk);
163     track->GetPxPyPz(momentum);
164     pt = TMath::Sqrt(momentum[0]*momentum[0]+momentum[1]*momentum[1]);
165
166     phitrk = TMath::ATan2(postrk[1]-posVtx[1],postrk[0]-posVtx[0]);
167     sigmavtxphi = TMath::Sqrt(sigmaVtx[0]*sigmaVtx[0]*
168                               TMath::Cos(phitrk)*TMath::Cos(phitrk)+
169                               sigmaVtx[1]*sigmaVtx[1]*
170                               TMath::Sin(phitrk)*TMath::Sin(phitrk)+
171                               covVtxXY*
172                               TMath::Cos(phitrk)*TMath::Sin(phitrk));
173     sigma = TMath::Sqrt(Sigmad0rphi(pt)*Sigmad0rphi(pt)+
174                         sigmavtxphi*sigmavtxphi);
175     maxd0rphi = fNSigma*sigma;
176     if(OptImpParCut==1) maxd0rphi *= 5.;
177
178     if(fDebug) printf("trk %d; lab %d; |d0| = %f;  cut = %f\n",i,track->GetLabel(),TMath::Abs(d0z0[0]),maxd0rphi);
179     if(TMath::Abs(d0z0[0]) > maxd0rphi) { 
180       if(fDebug) printf("     rejected\n");
181       delete track; continue; 
182     }
183
184     fTrkArray.AddLast(track);
185     nTrks++; 
186   }
187
188   delete initVertex;
189
190   return nTrks;
191
192 //----------------------------------------------------------------------------
193 Double_t AliVertexerTracks::Sigmad0rphi(Double_t pt) const {
194 //
195 // Impact parameter resolution in rphi [cm] vs pt [GeV/c]
196 //
197   Double_t a_meas  = 11.6;
198   Double_t a_scatt = 65.8;
199   Double_t b       = 1.878;
200
201   Double_t sigma = a_meas*a_meas+a_scatt*a_scatt/TMath::Power(pt,b);
202   sigma = 0.0001*TMath::Sqrt(sigma);
203
204   return sigma;
205 }
206 //----------------------------------------------------------------------------
207 AliVertex* AliVertexerTracks::VertexForSelectedTracks(TTree *trkTree) {
208 //
209 // Return vertex from tracks in trkTree
210 //
211
212   // get tracks and propagate them to initial vertex position
213   Int_t nTrks = PrepareTracks(*trkTree,1);
214   if(nTrks < fMinTracks) {
215     printf("TooFewTracks\n");
216     Double_t vtx[3]={0,0,0};
217     fVert.SetXYZ(vtx);
218     fVert.SetDispersion(999);
219     fVert.SetNContributors(-5);
220     return &fVert;
221   }
222  
223   // Set initial vertex position from ESD
224   if(fAlgo==1)  StrLinVertexFinderMinDist(1);
225   if(fAlgo==2)  StrLinVertexFinderMinDist(0);
226   if(fAlgo==3)  HelixVertexFinder();
227   if(fAlgo==4)  VertexFinder(1);
228   if(fAlgo==5)  VertexFinder(0);
229   return &fVert;
230 }
231 //----------------------------------------------------------------------------
232 AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESD *esdEvent)
233 {
234 //
235 // Primary vertex for current ESD event
236 // (Two iterations: 
237 //  1st with 5*fNSigma*sigma(pt) cut w.r.t. to initial vertex; 
238 //  2nd with fNSigma*sigma(pt) cut w.r.t. to vertex found in 1st iteration) 
239 //
240   fCurrentVertex = 0;
241
242   Int_t entr = (Int_t)esdEvent->GetNumberOfTracks();
243   TTree *trkTree = new TTree("TreeT","tracks");
244   AliESDtrack *esdTrack = 0;
245   trkTree->Branch("tracks","AliESDtrack",&esdTrack);
246
247   Bool_t   skipThis;
248   for(Int_t i=0; i<entr; i++) {
249     // check tracks to skip
250     skipThis = kFALSE;
251     for(Int_t j=0; j<fNTrksToSkip; j++) { 
252       if(i==fTrksToSkip[j]) {
253         if(fDebug) printf("skipping track: %d\n",i);
254         skipThis = kTRUE;
255       }
256     }
257     if(skipThis) continue;
258     AliESDtrack *et = esdEvent->GetTrack(i);
259     esdTrack = new AliESDtrack(*et);
260     if(!esdTrack->GetStatus()&AliESDtrack::kITSin) continue;
261     if(!esdTrack->GetStatus()&AliESDtrack::kITSrefit) continue;
262     Int_t nclus=esdTrack->GetNcls(0); // check number of clusters in ITS
263     if(nclus<fMinITSClusters) continue;
264     trkTree->Fill();
265   }
266   delete esdTrack;
267
268   // ITERATION 1
269   // propagate tracks to initVertex
270   // preselect them  (reject for |d0|>5*fNSigma*sigma w.r.t. initVertex)
271   Int_t nTrks;
272   nTrks = PrepareTracks(*trkTree,1);
273   if(fDebug) printf(" tracks prepared - step 1: %d\n",nTrks);
274   if(nTrks < fMinTracks) {
275     printf("TooFewTracks\n");
276     TooFewTracks(esdEvent);
277     if(fDebug) fCurrentVertex->PrintStatus();
278     return fCurrentVertex; 
279   }
280
281   // vertex finder
282   switch (fAlgo) {
283     case 1: StrLinVertexFinderMinDist(1); break;
284     case 2: StrLinVertexFinderMinDist(0); break;
285     case 3: HelixVertexFinder();          break;
286     case 4: VertexFinder(1);              break;
287     case 5: VertexFinder(0);              break;
288     default: printf("Wrong algorithm\n"); break;  
289   }
290   if(fDebug) printf(" vertex finding completed\n");
291
292   // vertex fitter
293   VertexFitter(kTRUE);
294   if(fDebug) printf(" vertex fit completed\n");
295
296   // ITERATION 2
297   // propagate tracks to best between initVertex and fCurrentVertex
298   // preselect tracks (reject for |d0|>fNSigma*sigma w.r.t. best 
299   //                   between initVertex and fCurrentVertex) 
300   nTrks = PrepareTracks(*trkTree,2);
301   delete trkTree;
302   if(fDebug) printf(" tracks prepared - step 2: %d\n",nTrks);
303   if(nTrks < fMinTracks) {
304     printf("TooFewTracks\n");
305     TooFewTracks(esdEvent);
306     if(fDebug) fCurrentVertex->PrintStatus();
307     return fCurrentVertex; 
308   }
309
310   // vertex finder
311   switch (fAlgo) {
312     case 1: StrLinVertexFinderMinDist(1); break;
313     case 2: StrLinVertexFinderMinDist(0); break;
314     case 3: HelixVertexFinder();          break;
315     case 4: VertexFinder(1);              break;
316     case 5: VertexFinder(0);              break;
317     default: printf("Wrong algorithm\n"); break;  
318   }
319   if(fDebug) printf(" vertex finding completed\n");
320
321   // fitter
322   VertexFitter(kTRUE);
323   if(fDebug) printf(" vertex fit completed\n");
324
325
326   fTrkArray.Clear();
327
328   // take true pos from ESD
329   Double_t tp[3];
330   esdEvent->GetVertex()->GetTruePos(tp);
331   fCurrentVertex->SetTruePos(tp);
332   if(fNominalSigma[0]>1.) {
333     fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
334   } else {
335     fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
336   }
337
338   if(fDebug) fCurrentVertex->PrintStatus();
339   if(fTrksToSkip) delete [] fTrksToSkip;
340
341   return fCurrentVertex;
342 }
343 //----------------------------------------------------------------------------
344 AliESDVertex* AliVertexerTracks::FindPrimaryVertexOld(const AliESD *esdEvent)
345 {
346 //
347 // Primary vertex for current ESD event (one iteration with |d0rphi|<3cm)
348 //
349   fCurrentVertex = 0;
350
351   Int_t entr = (Int_t)esdEvent->GetNumberOfTracks();
352   TTree *trkTree = new TTree("TreeT","tracks");
353   AliESDtrack *esdTrack = 0;
354   trkTree->Branch("tracks","AliESDtrack",&esdTrack);
355
356   for(Int_t i=0; i<entr; i++) {
357     AliESDtrack *et = esdEvent->GetTrack(i);
358     esdTrack = new AliESDtrack(*et);
359     if(!esdTrack->GetStatus()&AliESDtrack::kITSin) continue;
360     if(!esdTrack->GetStatus()&AliESDtrack::kITSrefit) continue;
361     Int_t nclus=esdTrack->GetNcls(0); // check number of clusters in ITS
362     if(nclus<fMinITSClusters) continue;
363
364     trkTree->Fill();
365   }
366   delete esdTrack;
367
368   // preselect tracks and propagate them to initial vertex position
369   Int_t nTrks = PrepareTracks(*trkTree,1);
370   delete trkTree;
371   if(fDebug) printf(" tracks prepared: %d\n",nTrks);
372   if(nTrks < fMinTracks) {
373     printf("TooFewTracks\n");
374     fCurrentVertex = new AliESDVertex(0.,0.,-1);
375     return fCurrentVertex; 
376   }
377
378   // VERTEX FINDER
379   switch (fAlgo) {
380     case 1: StrLinVertexFinderMinDist(1); break;
381     case 2: StrLinVertexFinderMinDist(0); break;
382     case 3: HelixVertexFinder();          break;
383     case 4: VertexFinder(1);              break;
384     case 5: VertexFinder(0);              break;
385     default: printf("Wrong algorithm\n"); break;  
386   }
387
388   // VERTEX FITTER
389   VertexFitter();
390   if(fDebug) printf(" vertex fit completed\n");
391
392
393   Double_t tp[3];
394   esdEvent->GetVertex()->GetTruePos(tp);
395   fCurrentVertex->SetTruePos(tp);
396   fCurrentVertex->SetTitle("VertexerTracks");
397
398   fTrkArray.Clear();
399   return fCurrentVertex;
400 }
401 //----------------------------------------------------------------------------
402 AliVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray) {
403 //
404 // Return vertex from array of tracks
405 //
406
407   // get tracks and propagate them to initial vertex position
408   Int_t nTrks = trkArray->GetEntriesFast();
409   if(nTrks < fMinTracks) {
410     printf("TooFewTracks\n");
411     Double_t vtx[3]={0,0,0};
412     fVert.SetXYZ(vtx);
413     fVert.SetDispersion(999);
414     fVert.SetNContributors(-5);
415     return &fVert;
416   }
417   TTree *trkTree = new TTree("TreeT","tracks");
418   AliESDtrack *esdTrack = 0;
419   trkTree->Branch("tracks","AliESDtrack",&esdTrack);
420   for(Int_t i=0; i<nTrks; i++){
421     esdTrack = (AliESDtrack*)trkArray->At(i);
422     trkTree->Fill();
423   }
424     
425   AliVertex *vtx =  VertexForSelectedTracks(trkTree);
426   delete trkTree;
427   return vtx;
428 }
429
430 //---------------------------------------------------------------------------
431 void AliVertexerTracks::VertexFinder(Int_t optUseWeights) {
432
433   // Get estimate of vertex position in (x,y) from tracks DCA
434  
435   Double_t initPos[3];
436   initPos[2] = 0.;
437   for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
438   Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
439   Double_t aver[3]={0.,0.,0.};
440   Double_t aversq[3]={0.,0.,0.};
441   Double_t sigmasq[3]={0.,0.,0.};
442   Double_t sigma=0;
443   Int_t ncombi = 0;
444   AliESDtrack *track1;
445   AliESDtrack *track2;
446   Double_t pos[3],dir[3]; 
447   Double_t alpha,mindist;
448   Double_t field=GetField();
449
450   for(Int_t i=0; i<nacc; i++){
451     track1 = (AliESDtrack*)fTrkArray.At(i);
452     alpha=track1->GetAlpha();
453     mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
454     track1->GetXYZAt(mindist,field,pos);
455     track1->GetPxPyPzAt(mindist,field,dir);
456     AliStrLine *line1 = new AliStrLine(pos,dir); 
457
458    //    AliStrLine *line1 = new AliStrLine();
459    //    track1->ApproximateHelixWithLine(mindist,field,line1);
460    
461     for(Int_t j=i+1; j<nacc; j++){
462       track2 = (AliESDtrack*)fTrkArray.At(j);
463       alpha=track2->GetAlpha();
464       mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
465       track2->GetXYZAt(mindist,field,pos);
466       track2->GetPxPyPzAt(mindist,field,dir);
467       AliStrLine *line2 = new AliStrLine(pos,dir); 
468     //      AliStrLine *line2 = new AliStrLine();
469     //  track2->ApproximateHelixWithLine(mindist,field,line2);
470       Double_t distCA=line2->GetDCA(line1);
471       if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
472         Double_t pnt1[3],pnt2[3],crosspoint[3];
473
474         if(optUseWeights<=0){
475           Int_t retcode = line2->Cross(line1,crosspoint);
476           if(retcode>=0){
477             ncombi++;
478             for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
479             for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
480           }
481         }
482         if(optUseWeights>0){
483           Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2);
484           if(retcode>=0){
485             Double_t alpha, cs, sn;
486             alpha=track1->GetAlpha();
487             cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);   
488             Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
489             alpha=track2->GetAlpha();
490             cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
491             Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
492             Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
493             Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
494             Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
495             Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
496             crosspoint[0]=wx1*pnt1[0] + wx2*pnt2[0]; 
497             crosspoint[1]=wy1*pnt1[1] + wy2*pnt2[1]; 
498             crosspoint[2]=wz1*pnt1[2] + wz2*pnt2[2];
499           
500             ncombi++;
501             for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
502             for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
503           }
504         }
505       }
506       delete line2;
507     }
508     delete line1;
509   }
510   if(ncombi>0){
511     for(Int_t jj=0;jj<3;jj++){
512       initPos[jj] = aver[jj]/ncombi;
513       aversq[jj]/=ncombi;
514       sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
515       sigma+=sigmasq[jj];
516     }
517     sigma=TMath::Sqrt(TMath::Abs(sigma));
518   }
519   else {
520     Warning("VertexFinder","Finder did not succed");
521     sigma=999;
522   }
523   fVert.SetXYZ(initPos);
524   fVert.SetDispersion(sigma);
525   fVert.SetNContributors(ncombi);
526 }
527 //---------------------------------------------------------------------------
528 void AliVertexerTracks::HelixVertexFinder() {
529
530   // Get estimate of vertex position in (x,y) from tracks DCA
531
532
533   Double_t initPos[3];
534   initPos[2] = 0.;
535   for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
536   Double_t field=GetField();
537
538   Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
539
540   Double_t aver[3]={0.,0.,0.};
541   Double_t averquad[3]={0.,0.,0.};
542   Double_t sigmaquad[3]={0.,0.,0.};
543   Double_t sigma=0;
544   Int_t ncombi = 0;
545   AliESDtrack *track1;
546   AliESDtrack *track2;
547   Double_t distCA;
548   Double_t x, par[5];
549   Double_t alpha, cs, sn;
550   Double_t crosspoint[3];
551   for(Int_t i=0; i<nacc; i++){
552     track1 = (AliESDtrack*)fTrkArray.At(i);
553     
554
555     for(Int_t j=i+1; j<nacc; j++){
556       track2 = (AliESDtrack*)fTrkArray.At(j);
557
558       distCA=track2->PropagateToDCA(track1,field);
559
560       if(fDCAcut<=0 ||(fDCAcut>0&&distCA<fDCAcut)){
561         track1->GetExternalParameters(x,par);
562         alpha=track1->GetAlpha();
563         cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
564         Double_t x1=x*cs - par[0]*sn;
565         Double_t y1=x*sn + par[0]*cs;
566         Double_t z1=par[1];
567         Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2(); 
568         track2->GetExternalParameters(x,par);
569         alpha=track2->GetAlpha();
570         cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
571         Double_t x2=x*cs - par[0]*sn;
572         Double_t y2=x*sn + par[0]*cs;
573         Double_t z2=par[1];
574         Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
575         Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
576         Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
577         Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
578         Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
579         crosspoint[0]=wx1*x1 + wx2*x2; 
580         crosspoint[1]=wy1*y1 + wy2*y2; 
581         crosspoint[2]=wz1*z1 + wz2*z2;
582
583         ncombi++;
584         for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
585         for(Int_t jj=0;jj<3;jj++)averquad[jj]+=(crosspoint[jj]*crosspoint[jj]);
586       }
587     }
588       
589   }
590   if(ncombi>0){
591     for(Int_t jj=0;jj<3;jj++){
592       initPos[jj] = aver[jj]/ncombi;
593       averquad[jj]/=ncombi;
594       sigmaquad[jj]=averquad[jj]-initPos[jj]*initPos[jj];
595       sigma+=sigmaquad[jj];
596     }
597     sigma=TMath::Sqrt(TMath::Abs(sigma));
598   }
599   else {
600     Warning("HelixVertexFinder","Finder did not succed");
601     sigma=999;
602   }
603   fVert.SetXYZ(initPos);
604   fVert.SetDispersion(sigma);
605   fVert.SetNContributors(ncombi);
606 }
607 //---------------------------------------------------------------------------
608 void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights){
609
610   // Calculate the point at minimum distance to prepared tracks 
611   
612   Double_t initPos[3];
613   initPos[2] = 0.;
614   Double_t sigma=0;
615   for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
616   const Int_t knacc = (Int_t)fTrkArray.GetEntriesFast();
617   Double_t field=GetField();
618
619   AliESDtrack *track1;
620   Double_t (*vectP0)[3]=new Double_t [knacc][3];
621   Double_t (*vectP1)[3]=new Double_t [knacc][3];
622   
623   Double_t sum[3][3];
624   Double_t dsum[3]={0,0,0};
625   for(Int_t i=0;i<3;i++)
626     for(Int_t j=0;j<3;j++)sum[i][j]=0;
627   for(Int_t i=0; i<knacc; i++){
628     track1 = (AliESDtrack*)fTrkArray.At(i);
629     Double_t alpha=track1->GetAlpha();
630     Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
631     Double_t pos[3],dir[3]; 
632     track1->GetXYZAt(mindist,field,pos);
633     track1->GetPxPyPzAt(mindist,field,dir);
634     AliStrLine *line1 = new AliStrLine(pos,dir); 
635     //    AliStrLine *line1 = new AliStrLine();
636     //    track1->ApproximateHelixWithLine(mindist,field,line1);
637
638     Double_t p0[3],cd[3];
639     line1->GetP0(p0);
640     line1->GetCd(cd);
641     Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
642     vectP0[i][0]=p0[0];
643     vectP0[i][1]=p0[1];
644     vectP0[i][2]=p0[2];
645     vectP1[i][0]=p1[0];
646     vectP1[i][1]=p1[1];
647     vectP1[i][2]=p1[2];
648     
649     Double_t matr[3][3];
650     Double_t dknow[3];
651     if(optUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
652     if(optUseWeights==1){
653       Double_t sigmasq[3];
654       sigmasq[0]=track1->GetSigmaY2();
655       sigmasq[1]=track1->GetSigmaY2();
656       sigmasq[2]=track1->GetSigmaZ2();
657       GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
658     }
659
660     for(Int_t iii=0;iii<3;iii++){
661       dsum[iii]+=dknow[iii]; 
662       for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
663     }
664     delete line1;
665   }
666   
667   Double_t vett[3][3];
668   Double_t det=GetDeterminant3X3(sum);
669   
670    if(det!=0){
671      for(Int_t zz=0;zz<3;zz++){
672        for(Int_t ww=0;ww<3;ww++){
673          for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
674        }
675        for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
676        initPos[zz]=GetDeterminant3X3(vett)/det;
677      }
678
679
680      for(Int_t i=0; i<knacc; i++){
681        Double_t p0[3]={0,0,0},p1[3]={0,0,0};
682        for(Int_t ii=0;ii<3;ii++){
683          p0[ii]=vectP0[i][ii];
684          p1[ii]=vectP1[i][ii];
685        }
686        sigma+=GetStrLinMinDist(p0,p1,initPos);
687      }
688
689      sigma=TMath::Sqrt(sigma);
690    }else{
691     Warning("StrLinVertexFinderMinDist","Finder did not succed");
692     sigma=999;
693   }
694   delete vectP0;
695   delete vectP1;
696   fVert.SetXYZ(initPos);
697   fVert.SetDispersion(sigma);
698   fVert.SetNContributors(knacc);
699 }
700 //_______________________________________________________________________
701 Double_t AliVertexerTracks::GetDeterminant3X3(Double_t matr[][3]){
702   //
703   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];
704  return det;
705 }
706 //____________________________________________________________________________
707 void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t (*m)[3],Double_t *d){
708
709   //
710   Double_t x12=p0[0]-p1[0];
711   Double_t y12=p0[1]-p1[1];
712   Double_t z12=p0[2]-p1[2];
713   Double_t kk=x12*x12+y12*y12+z12*z12;
714   m[0][0]=2-2/kk*x12*x12;
715   m[0][1]=-2/kk*x12*y12;
716   m[0][2]=-2/kk*x12*z12;
717   m[1][0]=-2/kk*x12*y12;
718   m[1][1]=2-2/kk*y12*y12;
719   m[1][2]=-2/kk*y12*z12;
720   m[2][0]=-2/kk*x12*z12;
721   m[2][1]=-2*y12*z12;
722   m[2][2]=2-2/kk*z12*z12;
723   d[0]=2*p0[0]-2/kk*p0[0]*x12*x12-2/kk*p0[2]*x12*z12-2/kk*p0[1]*x12*y12;
724   d[1]=2*p0[1]-2/kk*p0[1]*y12*y12-2/kk*p0[0]*x12*y12-2/kk*p0[2]*z12*y12;
725   d[2]=2*p0[2]-2/kk*p0[2]*z12*z12-2/kk*p0[0]*x12*z12-2/kk*p0[1]*z12*y12;
726
727 }
728 //____________________________________________________________________________
729 void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t *sigmasq,Double_t (*m)[3],Double_t *d){
730   //
731   Double_t x12=p1[0]-p0[0];
732   Double_t y12=p1[1]-p0[1];
733   Double_t z12=p1[2]-p0[2];
734
735   Double_t den= x12*x12*sigmasq[1]*sigmasq[2]+y12*y12*sigmasq[0]*sigmasq[2]+z12*z12*sigmasq[0]*sigmasq[1];
736
737   Double_t kk= 2*(x12*x12/sigmasq[0]+y12*y12/sigmasq[1]+z12*z12/sigmasq[2]);
738
739   Double_t cc[3];
740   cc[0]=-x12/sigmasq[0];
741   cc[1]=-y12/sigmasq[1];
742   cc[2]=-z12/sigmasq[2];
743
744   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;
745
746   Double_t ss= -p0[0]*cc[0]-p0[1]*cc[1]-p0[2]*cc[2];
747
748   Double_t aa[3];
749   aa[0]=x12*sigmasq[1]*sigmasq[2]/den;
750   aa[1]=y12*sigmasq[0]*sigmasq[2]/den;
751   aa[2]=z12*sigmasq[0]*sigmasq[1]/den;
752
753   m[0][0]=aa[0]*(aa[0]*kk+2*cc[0])+2*cc[0]*aa[0]+2/sigmasq[0];
754   m[0][1]=aa[1]*(aa[0]*kk+2*cc[0])+2*cc[1]*aa[0];
755   m[0][2]=aa[2]*(aa[0]*kk+2*cc[0])+2*cc[2]*aa[0];
756
757   m[1][0]=aa[0]*(aa[1]*kk+2*cc[1])+2*cc[0]*aa[1];
758   m[1][1]=aa[1]*(aa[1]*kk+2*cc[1])+2*cc[1]*aa[1]+2/sigmasq[1];
759   m[1][2]=aa[2]*(aa[1]*kk+2*cc[1])+2*cc[2]*aa[1];
760
761   m[2][0]=aa[0]*(aa[2]*kk+2*cc[2])+2*cc[0]*aa[2];
762   m[2][1]=aa[1]*(aa[2]*kk+2*cc[2])+2*cc[1]*aa[2];
763   m[2][2]=aa[2]*(aa[2]*kk+2*cc[2])+2*cc[2]*aa[2]+2/sigmasq[2];
764
765   d[0]=-ww*(aa[0]*kk+2*cc[0])-2*ss*aa[0]+2*p0[0]/sigmasq[0];
766   d[1]=-ww*(aa[1]*kk+2*cc[1])-2*ss*aa[1]+2*p0[1]/sigmasq[1];
767   d[2]=-ww*(aa[2]*kk+2*cc[2])-2*ss*aa[2]+2*p0[2]/sigmasq[2];
768
769   }
770 //_____________________________________________________________________________
771 Double_t AliVertexerTracks::GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t *x0){
772   //
773   Double_t x12=p0[0]-p1[0];
774   Double_t y12=p0[1]-p1[1];
775   Double_t z12=p0[2]-p1[2];
776   Double_t x10=p0[0]-x0[0];
777   Double_t y10=p0[1]-x0[1];
778   Double_t z10=p0[2]-x0[2];
779   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);
780 }
781 //---------------------------------------------------------------------------
782 void AliVertexerTracks::VertexFitter(Bool_t useNominalVtx) {
783 //
784 // The optimal estimate of the vertex position is given by a "weighted 
785 // average of tracks positions"
786 // Original method: CMS Note 97/0051
787 //
788   Double_t initPos[3];
789   fVert.GetXYZ(initPos);
790   if(fDebug) { 
791     printf(" VertexFitter(): start\n");
792     printf(" Number of tracks in array: %d\n",(Int_t)fTrkArray.GetEntriesFast());
793     printf(" Minimum # tracks required in fit: %d\n",fMinTracks);
794     printf("Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]);
795     if(useNominalVtx) printf(" This vertex will be used in fit: (%f+-%f,%f+-%f)\n",fNominalPos[0],fNominalSigma[0],fNominalPos[1],fNominalSigma[1]); 
796   }
797
798   Int_t i,j,k,step=0;
799   TMatrixD rv(3,1);
800   TMatrixD vV(3,3);
801   rv(0,0) = initPos[0];
802   rv(1,0) = initPos[1];
803   rv(2,0) = 0.;
804   Double_t xlStart,alpha;
805   Double_t rotAngle;
806   Double_t cosRot,sinRot;
807   Double_t cc[15];
808   Int_t nUsedTrks;
809   Double_t chi2,chi2i;
810   Int_t arrEntries = (Int_t)fTrkArray.GetEntries();
811   AliESDtrack *t = 0;
812   Int_t failed = 0;
813
814   // 2 steps:
815   // 1st - estimate of vtx using all tracks
816   // 2nd - estimate of global chi2
817   for(step=0; step<2; step++) {
818     if(fDebug) printf(" step = %d\n",step);
819     chi2 = 0.;
820     nUsedTrks = 0;
821
822     TMatrixD sumWiri(3,1);
823     TMatrixD sumWi(3,3);
824     for(i=0; i<3; i++) {
825       sumWiri(i,0) = 0.;
826       for(j=0; j<3; j++) sumWi(j,i) = 0.;
827     }
828
829     if(useNominalVtx) {
830       for(i=0; i<3; i++) {
831         sumWiri(i,0) += fNominalPos[i]/fNominalSigma[i]/fNominalSigma[i];
832         sumWi(i,i)   += 1./fNominalSigma[i]/fNominalSigma[i];
833       }
834     }
835
836
837     // loop on tracks  
838     for(k=0; k<arrEntries; k++) {
839       // get track from track array
840       t = (AliESDtrack*)fTrkArray.At(k);
841       alpha = t->GetAlpha();
842       xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
843       t->PropagateTo(xlStart,AliTracker::GetBz());   // to vtxSeed
844       rotAngle = alpha;
845       if(alpha<0.) rotAngle += 2.*TMath::Pi();
846       cosRot = TMath::Cos(rotAngle);
847       sinRot = TMath::Sin(rotAngle);
848       
849       // vector of track global coordinates
850       TMatrixD ri(3,1);
851       ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
852       ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
853       ri(2,0) = t->GetZ();
854
855       // matrix to go from global (x,y,z) to local (y,z);
856       TMatrixD qQi(2,3);
857       qQi(0,0) = -sinRot;
858       qQi(0,1) = cosRot;
859       qQi(0,2) = 0.;
860       qQi(1,0) = 0.;
861       qQi(1,1) = 0.;
862       qQi(1,2) = 1.;
863
864       // covariance matrix of local (y,z) - inverted
865       TMatrixD uUi(2,2);
866       t->GetExternalCovariance(cc);
867       uUi(0,0) = cc[0];
868       uUi(0,1) = cc[1];
869       uUi(1,0) = cc[1];
870       uUi(1,1) = cc[2];
871
872       // weights matrix: wWi = qQiT * uUiInv * qQi
873       if(uUi.Determinant() <= 0.) continue;
874       TMatrixD uUiInv(TMatrixD::kInverted,uUi);
875       TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
876       TMatrixD wWi(qQi,TMatrixD::kTransposeMult,uUiInvQi);
877
878       // track chi2
879       TMatrixD deltar = rv; deltar -= ri;
880       TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
881       chi2i = deltar(0,0)*wWideltar(0,0)+
882               deltar(1,0)*wWideltar(1,0)+
883               deltar(2,0)*wWideltar(2,0);
884
885
886       // add to total chi2
887       chi2 += chi2i;
888
889       TMatrixD wWiri(wWi,TMatrixD::kMult,ri); 
890
891       sumWiri += wWiri;
892       sumWi   += wWi;
893
894       nUsedTrks++;
895     } // end loop on tracks
896
897     if(nUsedTrks < fMinTracks) {
898       failed=1;
899       continue;
900     }
901
902     Double_t determinant = sumWi.Determinant();
903     //cerr<<" determinant: "<<determinant<<endl;
904     if(determinant < 100.)  { 
905       printf("det(V) = 0\n");       
906       failed=1;
907       continue;
908     }
909
910     // inverted of weights matrix
911     TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
912     vV = invsumWi;
913      
914     // position of primary vertex
915     rv.Mult(vV,sumWiri);
916
917   } // end loop on the 2 steps
918
919   delete t;
920
921   if(failed) { 
922     printf("TooFewTracks\n");
923     fCurrentVertex = new AliESDVertex(0.,0.,-1);
924     return; 
925   }
926
927   Double_t position[3];
928   position[0] = rv(0,0);
929   position[1] = rv(1,0);
930   position[2] = rv(2,0);
931   Double_t covmatrix[6];
932   covmatrix[0] = vV(0,0);
933   covmatrix[1] = vV(0,1);
934   covmatrix[2] = vV(1,1);
935   covmatrix[3] = vV(0,2);
936   covmatrix[4] = vV(1,2);
937   covmatrix[5] = vV(2,2);
938   
939   // store data in the vertex object
940   fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nUsedTrks);
941
942   if(fDebug) {
943     printf(" VertexFitter(): finish\n");
944     printf(" rv = ( %f , %f , %f )\n\n",rv(0,0),rv(1,0),rv(2,0));
945     fCurrentVertex->PrintStatus();
946   }
947
948   return;
949 }
950 //---------------------------------------------------------------------------
951 void AliVertexerTracks::TooFewTracks(const AliESD* esdEvent) {
952 //
953 // When the number of tracks is < fMinTracks
954 //
955
956   // deal with vertices not found
957   Double_t pos[3],err[3];
958   Int_t    ncontr=0;
959   pos[0] = fNominalPos[0];
960   err[0] = fNominalSigma[0];
961   pos[1] = fNominalPos[1];
962   err[1] = fNominalSigma[1];
963   pos[2] = esdEvent->GetVertex()->GetZv();
964   err[2] = esdEvent->GetVertex()->GetZRes();
965   if(err[0]>1. && esdEvent->GetVertex()->GetNContributors()<=0) 
966     ncontr = -1; // (x,y,z) = (0,0,0) 
967   if(err[0]>1. && esdEvent->GetVertex()->GetNContributors()>0) 
968     ncontr = -2; // (x,y,z) = (0,0,z_fromSPD) 
969   if(err[0]<1. && esdEvent->GetVertex()->GetNContributors()<=0) 
970     ncontr = -3; // (x,y,z) = (x_mean,y_mean,0)
971   if(err[0]<1. && esdEvent->GetVertex()->GetNContributors()>0) 
972     ncontr = -4; // (x,y,z) = (x_mean,y_mean,z_fromSPD)
973   fCurrentVertex = 0;
974   fCurrentVertex = new AliESDVertex(pos,err);
975   fCurrentVertex->SetNContributors(ncontr);
976
977   Double_t tp[3];
978   esdEvent->GetVertex()->GetTruePos(tp);
979   fCurrentVertex->SetTruePos(tp);
980   fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
981   if(ncontr==-1||ncontr==-2) 
982     fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
983
984   return;
985 }