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