ITS/AliITSVertexerTracks.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 tracks
  18 //    It accepts V2 and ESD tracks
  19 //
  20 // Origin: A.Dainese, Padova,
  21 //         andrea.dainese@pd.infn.it
  22 //         M.Masera,  Torino,
  23 //         massimo.masera@to.infn.it
  24 //-----------------------------------------------------------------
  25
  26 //---- standard headers ----
  27 #include <Riostream.h>
  28 //---- Root headers --------
  29 #include <TKey.h>
  30 #include <TFile.h>
  31 #include <TTree.h>
  32 #include <TMatrixD.h>
  33 //---- AliRoot headers -----
  34 #include "AliITSStrLine.h"
  35 #include "AliITStrackV2.h"
  36 #include "AliESDVertex.h"
  37 #include "AliITSVertexerTracks.h"
  38 #include "AliESD.h"
  39 #include "AliESDtrack.h"
  40
  41
  42 ClassImp(AliITSVertexerTracks)
  43
  44
  45 //----------------------------------------------------------------------------
  46 AliITSVertexerTracks::AliITSVertexerTracks():AliITSVertexer() {
  47 //
  48 // Default constructor
  49 //
  50   fInFile  = 0;
  51   fOutFile = 0;
  52   SetVtxStart();
  53   SetMinTracks();
  54   fTrksToSkip = 0;
  55   fNTrksToSkip = 0;
  56   fDCAcut=0;
  57 }
  58 //----------------------------------------------------------------------------
  59 AliITSVertexerTracks::AliITSVertexerTracks(TFile *inFile,TFile *outFile,
  60                                            Int_t fEv,Int_t lEv,
  61                                            Double_t xStart,Double_t yStart) {
  62 //
  63 // Standard constructor
  64 //
  65   fCurrentVertex = 0;
  66   fInFile  = inFile;
  67   fOutFile = outFile;
  68   SetFirstEvent(fEv);
  69   SetLastEvent(lEv);
  70   SetVtxStart(xStart,yStart);
  71   SetMinTracks();
  72   fTrksToSkip = 0;
  73   fNTrksToSkip = 0;
  74   fDCAcut=0;
  75   SetDebug();
  76 }
  77 //----------------------------------------------------------------------------
  78 AliITSVertexerTracks::AliITSVertexerTracks(TString fn,
  79                                            Double_t xStart,Double_t yStart)
  80                                           :AliITSVertexer(fn) {
  81 //
  82 // Alternative constructor
  83 //
  84   fInFile  = 0;
  85   fOutFile = 0;
  86   SetVtxStart(xStart,yStart);
  87   SetMinTracks();
  88   fTrksToSkip = 0;
  89   fNTrksToSkip = 0;
  90   fDCAcut=0;
  91 }
  92 //______________________________________________________________________
  93 AliITSVertexerTracks::AliITSVertexerTracks(const AliITSVertexerTracks &vtxr) : AliITSVertexer(vtxr) {
  94   // Copy constructor
  95   // Copies are not allowed. The method is protected to avoid misuse.
  96   Error("AliITSVertexerTracks","Copy constructor not allowed\n");
  97 }
  98
  99 //______________________________________________________________________
 100 AliITSVertexerTracks& AliITSVertexerTracks::operator=(const AliITSVertexerTracks& /* vtxr */){
 101   // Assignment operator
 102   // Assignment is not allowed. The method is protected to avoid misuse.
 103   Error("= operator","Assignment operator not allowed\n");
 104   return *this;
 105 }
 106
 107 //-----------------------------------------------------------------------------
 108 AliITSVertexerTracks::~AliITSVertexerTracks() {
 109   // Default Destructor
 110   // The objects poited by the following pointers are not owned
 111   // by this class and are not deleted
 112
 113   fCurrentVertex = 0;
 114   fInFile        = 0;
 115   fOutFile       = 0;
 116 }
 117 //-----------------------------------------------------------------------------
 118 Bool_t AliITSVertexerTracks::CheckField() const {
 119 //
 120 // Check if the conv. const. has been set
 121 //
 122   AliITStrackV2 t;
 123   Double_t cc    = t.GetConvConst();
 124   Double_t field = 100./0.299792458/cc;
 125
 126   if(field<0.1 || field>0.6) {
 127     printf("AliITSVertexerTracks::CheckField():\n ERROR: AliKalmanTrack::fConvConst not set\n Use AliKalmanTrack::SetConvConst() or AliITSVertexerTracks::SetField()\n");
 128     return kFALSE;
 129   }
 130   printf("AliITSVertexerTracks::CheckField():  Using B = %3.1f T\n",field);
 131   return kTRUE;
 132 }
 133 //---------------------------------------------------------------------------
 134 void AliITSVertexerTracks::ComputeMaxChi2PerTrack(Int_t nTracks) {
 135 //
 136 // Max. contr. to the chi2 has been tuned as a function of multiplicity
 137 //
 138   if(nTracks < 7) { fMaxChi2PerTrack = 1.e6;
 139   } else { fMaxChi2PerTrack = 100.; }
 140
 141   return;
 142 }
 143 //---------------------------------------------------------------------------
 144 void AliITSVertexerTracks::FindVertices() {
 145 //
 146 // Vertices for all events from fFirstEvent to fLastEvent
 147 //
 148
 149   // Check if the conv. const. has been set
 150   if(!CheckField()) return;
 151
 152   TDirectory *curdir = 0;
 153
 154   // loop over events
 155   for(Int_t ev=fFirstEvent; ev<=fLastEvent; ev++) {
 156     if(ev % 100 == 0 || fDebug) printf("--- Processing event %d of %d ---\n",ev,fLastEvent);
 157
 158     FindPrimaryVertexForCurrentEvent(ev);
 159
 160     if(!fCurrentVertex) {
 161       printf("AliITSVertexerTracks::FindVertixes(): no tracks tree for event %d\n",ev);
 162       continue;
 163     }
 164
 165     if(fDebug) fCurrentVertex->PrintStatus();
 166
 167     // write vertex to file
 168     TString vtxName = "Vertex_";
 169     vtxName += ev;
 170     fCurrentVertex->SetName(vtxName.Data());
 171     fCurrentVertex->SetTitle("VertexerTracks");
 172     //WriteCurrentVertex();
 173     curdir = gDirectory;
 174     fOutFile->cd();
 175     fCurrentVertex->Write();
 176     curdir->cd();
 177     fCurrentVertex = 0;
 178   } // loop over events
 179
 180   return;
 181 }
 182 //---------------------------------------------------------------------------
 183 void AliITSVertexerTracks::FindVerticesESD() {
 184 //
 185 // Vertices for all events from fFirstEvent to fLastEvent
 186 //
 187
 188  // Check if the conv. const. has been set
 189  if(!CheckField()) return;
 190
 191  TDirectory *curdir = 0;
 192
 193  fInFile->cd();
 194  TTree *esdTree = (TTree*)fInFile->Get("esdTree");
 195  if(!esdTree) {
 196      printf("AliITSVertexerTracks::FindVerticesESD(): no tree in file!\n");
 197    return;
 198  }
 199  Int_t nev = (Int_t)esdTree->GetEntries();
 200  Int_t ev;
 201  // loop on events in tree
 202  for(Int_t i=0; i<nev; i++) {
 203    AliESD *esdEvent = new AliESD;
 204    esdTree->SetBranchAddress("ESD",&esdEvent);
 205    if(!esdTree->GetEvent(i)) {
 206      printf("AliITSVertexerTracks::FindVerticesESD(): not an ESD!\n");
 207      delete esdEvent;
 208      return;
 209    }
 210    ev = (Int_t)esdEvent->GetEventNumber();
 211    if(ev<fFirstEvent || ev>fLastEvent) { delete esdEvent; continue; }
 212    if(ev % 100 == 0 || fDebug)
 213      printf("--- Processing event %d of %d ---\n",ev,fLastEvent-fFirstEvent);
 214
 215    FindPrimaryVertexForCurrentEvent(esdEvent);
 216
 217    if(!fCurrentVertex) {
 218      printf("AliITSVertexerTracks::FindVertixesESD():\n no vertex for event %d\n",ev);
 219      continue;
 220    }
 221
 222    if(fDebug) fCurrentVertex->PrintStatus();
 223
 224    // write vertex to file
 225    TString vtxName = "Vertex_";
 226    vtxName += ev;
 227    fCurrentVertex->SetName(vtxName.Data());
 228    fCurrentVertex->SetTitle("VertexerTracks");
 229    //WriteCurrentVertex();
 230    curdir = gDirectory;
 231    fOutFile->cd();
 232    fCurrentVertex->Write();
 233    curdir->cd();
 234    fCurrentVertex = 0;
 235    esdEvent = 0;
 236    delete esdEvent;
 237  } // end loop over events
 238
 239  return;
 240 }
 241 //----------------------------------------------------------------------------
 242 Int_t AliITSVertexerTracks::PrepareTracks(TTree &trkTree) {
 243   //
 244   // Propagate tracks to initial vertex position and store them in a TObjArray
 245   //
 246   Double_t maxd0rphi = 3.;
 247   Int_t    nTrks    = 0;
 248   Bool_t   skipThis;
 249   Double_t d0rphi;
 250   Int_t    nEntries = (Int_t)trkTree.GetEntries();
 251
 252   if(!fTrkArray.IsEmpty()) fTrkArray.Clear();
 253   fTrkArray.Expand(nEntries);
 254
 255   if(fDebug) {
 256     printf(" PrepareTracks()\n");
 257     trkTree.Print();
 258   }
 259   cout<<" entr tree its tracks = "<<nEntries<<endl;
 260   for(Int_t i=0; i<nEntries; i++) {
 261     // check tracks to skip
 262     skipThis = kFALSE;
 263     for(Int_t j=0; j<fNTrksToSkip; j++) {
 264       if(i==fTrksToSkip[j]) {
 265         if(fDebug) printf("skipping track: %d\n",i);
 266         skipThis = kTRUE;
 267       }
 268     }
 269     if(skipThis) continue;
 270
 271     AliITStrackV2 *itstrack = new AliITStrackV2;
 272     trkTree.SetBranchAddress("tracks",&itstrack);
 273     trkTree.GetEvent(i);
 274     d0rphi=Prepare(itstrack);
 275     if(d0rphi> maxd0rphi) { if(fDebug)cout<<"    !!!! d0rphi "<<d0rphi<<endl;continue; }
 276     fTrkArray.AddLast(itstrack);
 277
 278     nTrks++;
 279     if(fDebug)cout<<" :-) nTrks, d0rphi "<<nTrks<<"  "<<d0rphi<<endl;
 280
 281   }
 282   if(fTrksToSkip) delete [] fTrksToSkip;
 283
 284   return nTrks;
 285 }
 286 //----------------------------------------------------------------------------
 287 Int_t AliITSVertexerTracks::PrepareTracks(AliESD* esdEvent,Int_t nofCand, Int_t *trkPos) {
 288   //
 289   // Propagate tracks to initial vertex position and store them in a TObjArray
 290   //
 291   Int_t    nTrks    = 0;
 292   Double_t maxd0rphi = 3.;
 293   Double_t d0rphi;
 294
 295   if(!fTrkArray.IsEmpty()) fTrkArray.Clear();
 296   fTrkArray.Expand(100);
 297
 298   if(fDebug) {
 299     printf(" PrepareTracks()\n");
 300   }
 301   AliITStrackV2* itstrack;
 302
 303   for(Int_t i=0; i<nofCand;i++){
 304     AliESDtrack *esdTrack = (AliESDtrack*)esdEvent->GetTrack(trkPos[i]);
 305     UInt_t status=esdTrack->GetStatus();
 306     if ((status&AliESDtrack::kTPCin)==0)continue;
 307     if ((status&AliESDtrack::kITSin)==0)continue;
 308     if ((status&AliESDtrack::kITSrefit)==0) continue;
 309
 310     itstrack = new AliITStrackV2(*esdTrack);
 311     d0rphi=Prepare(itstrack);
 312     if(d0rphi> maxd0rphi) { if(fDebug)cout<<"    !!!! d0rphi "<<d0rphi<<endl;continue; }
 313     Int_t nclus=itstrack->GetNumberOfClusters();
 314
 315     if(nclus<6){delete itstrack;continue;}
 316     fTrkArray.AddLast(itstrack);
 317
 318     nTrks++;
 319     if(fDebug)cout<<" :-) nTrks, d0rphi "<<nTrks<<"  "<<d0rphi<<endl;
 320     //delete itstrack;
 321   }
 322
 323
 324   if(fTrksToSkip) delete [] fTrksToSkip;
 325   return nTrks;
 326 }
 327 //----------------------------------------------------------------------------
 328 Double_t AliITSVertexerTracks::Prepare(AliITStrackV2* itstrack){
 329
 330   //
 331   Double_t alpha,xlStart,d0rphi;
 332   // propagate track to vtxSeed
 333   alpha  = itstrack->GetAlpha();
 334   xlStart = fNominalPos[0]*TMath::Cos(alpha)+fNominalPos[1]*TMath::Sin(alpha);
 335   if(itstrack->GetX()>3.)itstrack->PropagateTo(3.,0.0023,65.19); // to beam pipe (0.8 mm of Be)
 336   itstrack->PropagateTo(xlStart,0.,0.);
 337   // select tracks with d0rphi < maxd0rphi
 338   d0rphi = TMath::Abs(itstrack->GetD(fNominalPos[0],fNominalPos[1]));
 339   return d0rphi;
 340
 341 }
 342 //----------------------------------------------------------------------------
 343 void AliITSVertexerTracks::PrintStatus() const {
 344 //
 345 // Print status
 346 //
 347   printf(" Initial position (%f,%f)\n",fNominalPos[0],fNominalPos[1]);
 348   printf(" Vertex position after vertex finder:\n");
 349   fSimpVert.Print();
 350   printf(" Number of tracks in array: %d\n",(Int_t)fTrkArray.GetEntriesFast());
 351   printf(" Minimum # tracks required in fit: %d\n",fMinTracks);
 352
 353   return;
 354 }
 355 //----------------------------------------------------------------------------
 356 AliESDVertex* AliITSVertexerTracks::FindPrimaryVertexForCurrentEvent(Int_t evnumb) {
 357 //
 358 // Vertex for current event
 359 //
 360   fCurrentVertex = 0;
 361
 362   // get tree with tracks from input file
 363   TString treeName = "TreeT_ITS_";
 364   treeName += evnumb;
 365   TTree *trkTree=(TTree*)fInFile->Get(treeName.Data());
 366   if(!trkTree) return fCurrentVertex;
 367
 368
 369   // get tracks and propagate them to initial vertex position
 370   Int_t nTrks = PrepareTracks(*trkTree);
 371   delete trkTree;
 372   if(fDebug) printf(" tracks prepared: %d\n",nTrks);
 373   if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
 374
 375   // VERTEX FINDER
 376   VertexFinder();
 377
 378   // VERTEX FITTER
 379   ComputeMaxChi2PerTrack(nTrks);
 380   VertexFitter();
 381   if(fDebug) printf(" vertex fit completed\n");
 382
 383   return fCurrentVertex;
 384 }
 385 //----------------------------------------------------------------------------
 386 AliESDVertex* AliITSVertexerTracks::FindPrimaryVertexForCurrentEvent(AliESD *esdEvent)
 387 {
 388 //
 389 // Vertex for current ESD event
 390 //
 391   fCurrentVertex = 0;
 392   Double_t vtx[3],cvtx[6];
 393
 394   // put tracks reco in ITS in a tree
 395   Int_t entr = (Int_t)esdEvent->GetNumberOfTracks();
 396   TTree *trkTree = new TTree("TreeT_ITS","its tracks");
 397   AliITStrackV2 *itstrack = 0;
 398   trkTree->Branch("tracks","AliITStrackV2",&itstrack,entr,0);
 399
 400   for(Int_t i=0; i<entr; i++) {
 401     AliESDtrack *esdTrack = (AliESDtrack*)esdEvent->GetTrack(i);
 402     if(!esdTrack->GetStatus()&AliESDtrack::kITSin)
 403       { delete esdTrack; continue; }
 404     try {
 405       itstrack = new AliITStrackV2(*esdTrack);
 406     }
 407     catch (const Char_t *msg) {
 408         Warning("FindPrimaryVertexForCurrentEvent",msg);
 409         delete esdTrack;
 410         continue;
 411     }
 412
 413     trkTree->Fill();
 414     itstrack = 0;
 415     delete esdTrack;
 416   }
 417   delete itstrack;
 418
 419   // preselect tracks and propagate them to initial vertex position
 420   Int_t nTrks = PrepareTracks(*trkTree);
 421   delete trkTree;
 422   if(fDebug) printf(" tracks prepared: %d\n",nTrks);
 423   if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
 424
 425   // Set initial vertex position from ESD
 426   esdEvent->GetVertex()->GetXYZ(vtx);
 427   SetVtxStart(vtx[0],vtx[1]);
 428   // VERTEX FINDER
 429   VertexFinder();
 430
 431   // VERTEX FITTER
 432   ComputeMaxChi2PerTrack(nTrks);
 433   VertexFitter();
 434   if(fDebug) printf(" vertex fit completed\n");
 435
 436   // store vertex information in ESD
 437   fCurrentVertex->GetXYZ(vtx);
 438   fCurrentVertex->GetCovMatrix(cvtx);
 439
 440   Double_t tp[3];
 441   esdEvent->GetVertex()->GetTruePos(tp);
 442   fCurrentVertex->SetTruePos(tp);
 443
 444   esdEvent->SetVertex(fCurrentVertex);
 445
 446   cout<<"Vertex: "<<vtx[0]<<", "<<vtx[1]<<", "<<vtx[2]<<endl;
 447   return fCurrentVertex;
 448 }
 449 //----------------------------------------------------------------------------
 450 AliITSSimpleVertex* AliITSVertexerTracks::VertexForSelectedTracks(AliESD *esdEvent,Int_t nofCand, Int_t *trkPos,  Int_t opt){
 451
 452   //
 453   // Computes the vertex for selected tracks
 454   // trkPos=vector with track positions in ESD
 455   // values of opt -> see AliITSVertexerTracks.h
 456   //
 457   Double_t vtx[3]={0,0,0};
 458
 459   Int_t nTrks = PrepareTracks(esdEvent,nofCand, trkPos);
 460   //delete trkTree;//  :-))
 461   if(fDebug) printf(" tracks prepared: %d\n",nTrks);
 462   if(nTrks < fMinTracks) {
 463     fSimpVert.SetXYZ(vtx);
 464     fSimpVert.SetDispersion(999);
 465     fSimpVert.SetNContributors(-5);
 466     return &fSimpVert;
 467   }
 468
 469   // Set initial vertex position from ESD
 470   esdEvent->GetVertex()->GetXYZ(vtx);
 471   SetVtxStart(vtx[0],vtx[1]);
 472   if(opt==1)  StrLinVertexFinderMinDist(1);
 473   if(opt==2)  StrLinVertexFinderMinDist(0);
 474   if(opt==3)  HelixVertexFinder();
 475   if(opt==4)  VertexFinder(1);
 476   if(opt==5)  VertexFinder(0);
 477   return &fSimpVert;
 478 }
 479
 480 //---------------------------------------------------------------------------
 481 void AliITSVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) {
 482 //
 483 // Mark the tracks not ot be used in the vertex finding
 484 //
 485   fNTrksToSkip = n;
 486   fTrksToSkip = new Int_t[n];
 487   for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
 488   return;
 489 }
 490 //---------------------------------------------------------------------------
 491 void AliITSVertexerTracks::TooFewTracks() {
 492 //
 493 // When the number of tracks is < fMinTracks the vertex is set to (0,0,0)
 494 // and the number of tracks to -1
 495 //
 496   fCurrentVertex = new AliESDVertex(0.,0.,-1);
 497   return;
 498 }
 499 //---------------------------------------------------------------------------
 500 void AliITSVertexerTracks::VertexFinder(Int_t OptUseWeights) {
 501
 502   // Get estimate of vertex position in (x,y) from tracks DCA
 503   // Then this estimate is stored to the data member fSimpVert
 504   // (previous values are overwritten)
 505
 506
 507   Double_t initPos[3];
 508   initPos[2] = 0.;
 509   for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
 510   Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
 511   Double_t aver[3]={0.,0.,0.};
 512   Double_t aversq[3]={0.,0.,0.};
 513   Double_t sigmasq[3]={0.,0.,0.};
 514   Double_t sigma=0;
 515   Int_t ncombi = 0;
 516   AliITStrackV2 *track1;
 517   AliITStrackV2 *track2;
 518   Double_t alpha,mindist;
 519
 520   for(Int_t i=0; i<nacc; i++){
 521     track1 = (AliITStrackV2*)fTrkArray.At(i);
 522     alpha=track1->GetAlpha();
 523     mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
 524     AliITSStrLine *line1 = new AliITSStrLine();
 525     track1->ApproximateHelixWithLine(mindist,line1);
 526
 527     if(fDebug>5){
 528       Double_t xv,par[5];
 529       track1->GetExternalParameters(xv,par);
 530       cout<<"Track in position "<<i<<" xr= "<<xv<<endl;
 531       for(Int_t ii=0;ii<5;ii++)cout<<par[ii]<<" ";
 532       cout<<endl;
 533     }
 534
 535     for(Int_t j=i+1; j<nacc; j++){
 536       track2 = (AliITStrackV2*)fTrkArray.At(j);
 537       alpha=track2->GetAlpha();
 538       mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
 539       AliITSStrLine *line2 = new AliITSStrLine();
 540       track2->ApproximateHelixWithLine(mindist,line2);
 541       Double_t distCA=line2->GetDCA(line1);
 542       if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
 543         Double_t pnt1[3],pnt2[3],crosspoint[3];
 544
 545         if(OptUseWeights<=0){
 546           Int_t retcode = line2->Cross(line1,crosspoint);
 547           if(retcode<0){
 548             if(fDebug>10)cout<<" i= "<<i<<",   j= "<<j<<endl;
 549             if(fDebug>10)cout<<"bad intersection\n";
 550             line1->PrintStatus();
 551             line2->PrintStatus();
 552           }
 553           else {
 554             ncombi++;
 555             for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
 556             for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
 557             if(fDebug>10)cout<<" i= "<<i<<",   j= "<<j<<endl;
 558             if(fDebug>10)cout<<"\n Cross point: ";
 559             if(fDebug>10)cout<<crosspoint[0]<<" "<<crosspoint[1]<<" "<<crosspoint[2]<<endl;
 560           }
 561         }
 562         if(OptUseWeights>0){
 563           Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2);
 564           if(retcode>=0){
 565             Double_t alpha, cs, sn;
 566             alpha=track1->GetAlpha();
 567             cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
 568             Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
 569             alpha=track2->GetAlpha();
 570             cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
 571             Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
 572             Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
 573             Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
 574             Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
 575             Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
 576             crosspoint[0]=wx1*pnt1[0] + wx2*pnt2[0];
 577             crosspoint[1]=wy1*pnt1[1] + wy2*pnt2[1];
 578             crosspoint[2]=wz1*pnt1[2] + wz2*pnt2[2];
 579
 580             ncombi++;
 581             for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
 582             for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
 583           }
 584         }
 585       }
 586       delete line2;
 587     }
 588     delete line1;
 589   }
 590   if(ncombi>0){
 591     for(Int_t jj=0;jj<3;jj++){
 592       initPos[jj] = aver[jj]/ncombi;
 593       aversq[jj]/=ncombi;
 594       sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
 595       sigma+=sigmasq[jj];
 596     }
 597     sigma=TMath::Sqrt(TMath::Abs(sigma));
 598   }
 599   else {
 600     Warning("VertexFinder","Finder did not succed");
 601     sigma=999;
 602   }
 603   fSimpVert.SetXYZ(initPos);
 604   fSimpVert.SetDispersion(sigma);
 605   fSimpVert.SetNContributors(ncombi);
 606 }
 607 //---------------------------------------------------------------------------
 608 void AliITSVertexerTracks::HelixVertexFinder() {
 609
 610   // Get estimate of vertex position in (x,y) from tracks DCA
 611   // Then this estimate is stored to the data member fSimpVert
 612   // (previous values are overwritten)
 613
 614
 615   Double_t initPos[3];
 616   initPos[2] = 0.;
 617   for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
 618
 619   Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
 620
 621   Double_t aver[3]={0.,0.,0.};
 622   Double_t averquad[3]={0.,0.,0.};
 623   Double_t sigmaquad[3]={0.,0.,0.};
 624   Double_t sigma=0;
 625   Int_t ncombi = 0;
 626   AliITStrackV2 *track1;
 627   AliITStrackV2 *track2;
 628   Double_t distCA;
 629   Double_t x, par[5];
 630   Double_t alpha, cs, sn;
 631   Double_t crosspoint[3];
 632   for(Int_t i=0; i<nacc; i++){
 633     track1 = (AliITStrackV2*)fTrkArray.At(i);
 634
 635
 636     for(Int_t j=i+1; j<nacc; j++){
 637       track2 = (AliITStrackV2*)fTrkArray.At(j);
 638
 639
 640       distCA=track2->PropagateToDCA(track1);
 641
 642       if(fDCAcut<=0 ||(fDCAcut>0&&distCA<fDCAcut)){
 643         track1->GetExternalParameters(x,par);
 644         alpha=track1->GetAlpha();
 645         cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
 646         Double_t x1=x*cs - par[0]*sn;
 647         Double_t y1=x*sn + par[0]*cs;
 648         Double_t z1=par[1];
 649         Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
 650
 651         track2->GetExternalParameters(x,par);
 652         alpha=track2->GetAlpha();
 653         cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
 654         Double_t x2=x*cs - par[0]*sn;
 655         Double_t y2=x*sn + par[0]*cs;
 656         Double_t z2=par[1];
 657         Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
 658         //      printf("Track %d pos=(%f,%f,%f) - dca=%f\n",i,x1,y1,z1,distCA);
 659         //printf("Track %d pos=(%f,%f,%f)\n",j,x2,y2,z2);
 660
 661         Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
 662         Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
 663         Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
 664         Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
 665         crosspoint[0]=wx1*x1 + wx2*x2;
 666         crosspoint[1]=wy1*y1 + wy2*y2;
 667         crosspoint[2]=wz1*z1 + wz2*z2;
 668
 669         ncombi++;
 670         for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
 671         for(Int_t jj=0;jj<3;jj++)averquad[jj]+=(crosspoint[jj]*crosspoint[jj]);
 672       }
 673     }
 674
 675   }
 676   if(ncombi>0){
 677     for(Int_t jj=0;jj<3;jj++){
 678       initPos[jj] = aver[jj]/ncombi;
 679       averquad[jj]/=ncombi;
 680       sigmaquad[jj]=averquad[jj]-initPos[jj]*initPos[jj];
 681       sigma+=sigmaquad[jj];
 682     }
 683     sigma=TMath::Sqrt(TMath::Abs(sigma));
 684   }
 685   else {
 686     Warning("VertexFinder","Finder did not succed");
 687     sigma=999;
 688   }
 689   fSimpVert.SetXYZ(initPos);
 690   fSimpVert.SetDispersion(sigma);
 691   fSimpVert.SetNContributors(ncombi);
 692 }
 693 //---------------------------------------------------------------------------
 694 void AliITSVertexerTracks::StrLinVertexFinderMinDist(Int_t OptUseWeights){
 695
 696   // Calculate the point at minimum distance to prepared tracks
 697   // Then this estimate is stored to the data member fSimpVert
 698   // (previous values are overwritten)
 699
 700   Double_t initPos[3];
 701   initPos[2] = 0.;
 702   Double_t sigma=0;
 703   for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
 704   const Int_t knacc = (Int_t)fTrkArray.GetEntriesFast();
 705
 706   AliITStrackV2 *track1;
 707   Double_t (*vectP0)[3]=new Double_t [knacc][3];
 708   Double_t (*vectP1)[3]=new Double_t [knacc][3];
 709
 710   Double_t sum[3][3];
 711   Double_t dsum[3]={0,0,0};
 712   for(Int_t i=0;i<3;i++)
 713     for(Int_t j=0;j<3;j++)sum[i][j]=0;
 714   for(Int_t i=0; i<knacc; i++){
 715     track1 = (AliITStrackV2*)fTrkArray.At(i);
 716     Double_t alpha=track1->GetAlpha();
 717     Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
 718     AliITSStrLine *line1 = new AliITSStrLine();
 719     track1->ApproximateHelixWithLine(mindist,line1);
 720
 721     Double_t p0[3],cd[3];
 722     line1->GetP0(p0);
 723     line1->GetCd(cd);
 724     Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
 725     vectP0[i][0]=p0[0];
 726     vectP0[i][1]=p0[1];
 727     vectP0[i][2]=p0[2];
 728     vectP1[i][0]=p1[0];
 729     vectP1[i][1]=p1[1];
 730     vectP1[i][2]=p1[2];
 731
 732     Double_t matr[3][3];
 733     Double_t dknow[3];
 734     if(OptUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
 735     if(OptUseWeights==1){
 736       Double_t sigmasq[3];
 737       sigmasq[0]=track1->GetSigmaY2();
 738       sigmasq[1]=track1->GetSigmaY2();
 739       sigmasq[2]=track1->GetSigmaZ2();
 740       GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
 741     }
 742
 743     for(Int_t iii=0;iii<3;iii++){
 744       dsum[iii]+=dknow[iii];
 745       for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
 746     }
 747     delete line1;
 748   }
 749
 750   Double_t vett[3][3];
 751   Double_t det=GetDeterminant3X3(sum);
 752
 753    if(det!=0){
 754      for(Int_t zz=0;zz<3;zz++){
 755        for(Int_t ww=0;ww<3;ww++){
 756          for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
 757        }
 758        for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
 759        initPos[zz]=GetDeterminant3X3(vett)/det;
 760      }
 761
 762
 763      for(Int_t i=0; i<knacc; i++){
 764        Double_t p0[3]={0,0,0},p1[3]={0,0,0};
 765        for(Int_t ii=0;ii<3;ii++){
 766          p0[ii]=vectP0[i][ii];
 767          p1[ii]=vectP1[i][ii];
 768        }
 769        sigma+=GetStrLinMinDist(p0,p1,initPos);
 770      }
 771
 772      sigma=TMath::Sqrt(sigma);
 773    }else{
 774     Warning("VertexFinder","Finder did not succed");
 775     sigma=999;
 776   }
 777   delete vectP0;
 778   delete vectP1;
 779   fSimpVert.SetXYZ(initPos);
 780   fSimpVert.SetDispersion(sigma);
 781   fSimpVert.SetNContributors(knacc);
 782 }
 783 //_______________________________________________________________________
 784 Double_t AliITSVertexerTracks::GetDeterminant3X3(Double_t matr[][3]){
 785   //
 786   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];
 787  return det;
 788 }
 789 //____________________________________________________________________________
 790 void AliITSVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t m[][3],Double_t *d){
 791
 792   //
 793   Double_t x12=p0[0]-p1[0];
 794   Double_t y12=p0[1]-p1[1];
 795   Double_t z12=p0[2]-p1[2];
 796   Double_t kk=x12*x12+y12*y12+z12*z12;
 797   m[0][0]=2-2/kk*x12*x12;
 798   m[0][1]=-2/kk*x12*y12;
 799   m[0][2]=-2/kk*x12*z12;
 800   m[1][0]=-2/kk*x12*y12;
 801   m[1][1]=2-2/kk*y12*y12;
 802   m[1][2]=-2/kk*y12*z12;
 803   m[2][0]=-2/kk*x12*z12;
 804   m[2][1]=-2*y12*z12;
 805   m[2][2]=2-2/kk*z12*z12;
 806   d[0]=2*p0[0]-2/kk*p0[0]*x12*x12-2/kk*p0[2]*x12*z12-2/kk*p0[1]*x12*y12;
 807   d[1]=2*p0[1]-2/kk*p0[1]*y12*y12-2/kk*p0[0]*x12*y12-2/kk*p0[2]*z12*y12;
 808   d[2]=2*p0[2]-2/kk*p0[2]*z12*z12-2/kk*p0[0]*x12*z12-2/kk*p0[1]*z12*y12;
 809
 810 }
 811 //____________________________________________________________________________
 812 void AliITSVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t *sigmasq,Double_t m[][3],Double_t *d){
 813   //
 814   Double_t x12=p1[0]-p0[0];
 815   Double_t y12=p1[1]-p0[1];
 816   Double_t z12=p1[2]-p0[2];
 817
 818   Double_t den= x12*x12*sigmasq[1]*sigmasq[2]+y12*y12*sigmasq[0]*sigmasq[2]+z12*z12*sigmasq[0]*sigmasq[1];
 819
 820   Double_t kk= 2*(x12*x12/sigmasq[0]+y12*y12/sigmasq[1]+z12*z12/sigmasq[2]);
 821
 822   Double_t cc[3];
 823   cc[0]=-x12/sigmasq[0];
 824   cc[1]=-y12/sigmasq[1];
 825   cc[2]=-z12/sigmasq[2];
 826
 827   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;
 828
 829   Double_t ss= -p0[0]*cc[0]-p0[1]*cc[1]-p0[2]*cc[2];
 830
 831   Double_t aa[3];
 832   aa[0]=x12*sigmasq[1]*sigmasq[2]/den;
 833   aa[1]=y12*sigmasq[0]*sigmasq[2]/den;
 834   aa[2]=z12*sigmasq[0]*sigmasq[1]/den;
 835
 836   m[0][0]=aa[0]*(aa[0]*kk+2*cc[0])+2*cc[0]*aa[0]+2/sigmasq[0];
 837   m[0][1]=aa[1]*(aa[0]*kk+2*cc[0])+2*cc[1]*aa[0];
 838   m[0][2]=aa[2]*(aa[0]*kk+2*cc[0])+2*cc[2]*aa[0];
 839
 840   m[1][0]=aa[0]*(aa[1]*kk+2*cc[1])+2*cc[0]*aa[1];
 841   m[1][1]=aa[1]*(aa[1]*kk+2*cc[1])+2*cc[1]*aa[1]+2/sigmasq[1];
 842   m[1][2]=aa[2]*(aa[1]*kk+2*cc[1])+2*cc[2]*aa[1];
 843
 844   m[2][0]=aa[0]*(aa[2]*kk+2*cc[2])+2*cc[0]*aa[2];
 845   m[2][1]=aa[1]*(aa[2]*kk+2*cc[2])+2*cc[1]*aa[2];
 846   m[2][2]=aa[2]*(aa[2]*kk+2*cc[2])+2*cc[2]*aa[2]+2/sigmasq[2];
 847
 848   d[0]=-ww*(aa[0]*kk+2*cc[0])-2*ss*aa[0]+2*p0[0]/sigmasq[0];
 849   d[1]=-ww*(aa[1]*kk+2*cc[1])-2*ss*aa[1]+2*p0[1]/sigmasq[1];
 850   d[2]=-ww*(aa[2]*kk+2*cc[2])-2*ss*aa[2]+2*p0[2]/sigmasq[2];
 851
 852 }
 853 //_____________________________________________________________________________
 854 Double_t AliITSVertexerTracks::GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t *x0){
 855   //
 856   Double_t x12=p0[0]-p1[0];
 857   Double_t y12=p0[1]-p1[1];
 858   Double_t z12=p0[2]-p1[2];
 859   Double_t x10=p0[0]-x0[0];
 860   Double_t y10=p0[1]-x0[1];
 861   Double_t z10=p0[2]-x0[2];
 862   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);
 863
 864 }
 865 //---------------------------------------------------------------------------
 866 void AliITSVertexerTracks::VertexFitter() {
 867 //
 868 // The optimal estimate of the vertex position is given by a "weighted
 869 // average of tracks positions"
 870 // Original method: CMS Note 97/0051
 871 //
 872   if(fDebug) {
 873     printf(" VertexFitter(): start\n");
 874     PrintStatus();
 875   }
 876
 877
 878   Int_t i,j,k,step=0;
 879   TMatrixD rv(3,1);
 880   TMatrixD vV(3,3);
 881   Double_t initPos[3];
 882   fSimpVert.GetXYZ(initPos);
 883   rv(0,0) = initPos[0];
 884   rv(1,0) = initPos[1];
 885   rv(2,0) = 0.;
 886   Double_t xlStart,alpha;
 887   Double_t rotAngle;
 888   Double_t cosRot,sinRot;
 889   Double_t cc[15];
 890   Int_t nUsedTrks;
 891   Double_t chi2,chi2i;
 892   Int_t arrEntries = (Int_t)fTrkArray.GetEntries();
 893   AliITStrackV2 *t = 0;
 894   Int_t failed = 0;
 895
 896   Int_t *skipTrack = new Int_t[arrEntries];
 897   for(i=0; i<arrEntries; i++) skipTrack[i]=0;
 898
 899   // 3 steps:
 900   // 1st - first estimate of vtx using all tracks
 901   // 2nd - apply cut on chi2 max per track
 902   // 3rd - estimate of global chi2
 903   for(step=0; step<3; step++) {
 904     if(fDebug) printf(" step = %d\n",step);
 905     chi2 = 0.;
 906     nUsedTrks = 0;
 907
 908     TMatrixD sumWiri(3,1);
 909     TMatrixD sumWi(3,3);
 910     for(i=0; i<3; i++) {
 911       sumWiri(i,0) = 0.;
 912       for(j=0; j<3; j++) sumWi(j,i) = 0.;
 913     }
 914
 915     // loop on tracks
 916     for(k=0; k<arrEntries; k++) {
 917       if(skipTrack[k]) continue;
 918       // get track from track array
 919       t = (AliITStrackV2*)fTrkArray.At(k);
 920       alpha = t->GetAlpha();
 921       xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
 922       t->PropagateTo(xlStart,0.,0.);   // to vtxSeed
 923       rotAngle = alpha;
 924       if(alpha<0.) rotAngle += 2.*TMath::Pi();
 925       cosRot = TMath::Cos(rotAngle);
 926       sinRot = TMath::Sin(rotAngle);
 927
 928       // vector of track global coordinates
 929       TMatrixD ri(3,1);
 930       ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
 931       ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
 932       ri(2,0) = t->GetZ();
 933
 934       // matrix to go from global (x,y,z) to local (y,z);
 935       TMatrixD qQi(2,3);
 936       qQi(0,0) = -sinRot;
 937       qQi(0,1) = cosRot;
 938       qQi(0,2) = 0.;
 939       qQi(1,0) = 0.;
 940       qQi(1,1) = 0.;
 941       qQi(1,2) = 1.;
 942
 943       // covariance matrix of local (y,z) - inverted
 944       TMatrixD uUi(2,2);
 945       t->GetExternalCovariance(cc);
 946       uUi(0,0) = cc[0];
 947       uUi(0,1) = cc[1];
 948       uUi(1,0) = cc[1];
 949       uUi(1,1) = cc[2];
 950
 951       // weights matrix: wWi = qQiT * uUiInv * qQi
 952       if(uUi.Determinant() <= 0.) continue;
 953       TMatrixD uUiInv(TMatrixD::kInverted,uUi);
 954       TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
 955       TMatrixD wWi(qQi,TMatrixD::kTransposeMult,uUiInvQi);
 956
 957       // track chi2
 958       TMatrixD deltar = rv; deltar -= ri;
 959       TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
 960       chi2i = deltar(0,0)*wWideltar(0,0)+
 961               deltar(1,0)*wWideltar(1,0)+
 962               deltar(2,0)*wWideltar(2,0);
 963
 964
 965       if(step==1 && chi2i > fMaxChi2PerTrack) {
 966         skipTrack[k] = 1;
 967         continue;
 968       }
 969
 970       // add to total chi2
 971       chi2 += chi2i;
 972
 973       TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
 974
 975       sumWiri += wWiri;
 976       sumWi   += wWi;
 977
 978       nUsedTrks++;
 979     } // end loop on tracks
 980
 981     if(nUsedTrks < fMinTracks) {
 982       failed=1;
 983       continue;
 984     }
 985
 986     Double_t determinant = sumWi.Determinant();
 987     //cerr<<" determinant: "<<determinant<<endl;
 988     if(determinant < 100.)  {
 989       printf("det(V) = 0\n");
 990       failed=1;
 991       continue;
 992     }
 993
 994     // inverted of weights matrix
 995     TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
 996     vV = invsumWi;
 997
 998     // position of primary vertex
 999     rv.Mult(vV,sumWiri);
1000
1001   } // end loop on the 3 steps
1002
1003   delete [] skipTrack;
1004   delete t;
1005
1006   if(failed) {
1007     TooFewTracks();
1008     return;
1009   }
1010
1011   Double_t position[3];
1012   position[0] = rv(0,0);
1013   position[1] = rv(1,0);
1014   position[2] = rv(2,0);
1015   Double_t covmatrix[6];
1016   covmatrix[0] = vV(0,0);
1017   covmatrix[1] = vV(0,1);
1018   covmatrix[2] = vV(1,1);
1019   covmatrix[3] = vV(0,2);
1020   covmatrix[4] = vV(1,2);
1021   covmatrix[5] = vV(2,2);
1022
1023   // store data in the vertex object
1024   fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nUsedTrks);
1025
1026   if(fDebug) {
1027     printf(" VertexFitter(): finish\n");
1028     printf(" rv = ( %f , %f , %f )\n\n",rv(0,0),rv(1,0),rv(2,0));
1029     fCurrentVertex->PrintStatus();
1030   }
1031
1032   return;
1033 }
1034 //----------------------------------------------------------------------------
1035 AliESDVertex *AliITSVertexerTracks::VertexOnTheFly(TTree &trkTree) {
1036 //
1037 // Return vertex from tracks in trkTree
1038 //
1039   if(fCurrentVertex) fCurrentVertex = 0;
1040
1041   // get tracks and propagate them to initial vertex position
1042   Int_t nTrks = PrepareTracks(*(&trkTree));
1043   if(fDebug) printf(" tracks prepared: %d\n",nTrks);
1044   if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
1045
1046   // VERTEX FINDER
1047   VertexFinder();
1048
1049   // VERTEX FITTER
1050   ComputeMaxChi2PerTrack(nTrks);
1051   VertexFitter();
1052   if(fDebug) printf(" vertex fit completed\n");
1053
1054   return fCurrentVertex;
1055 }
1056 //----------------------------------------------------------------------------
1057
1058
1059