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