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