New vertex finders optimized for 3 prongs secondary vertices (AliITSVertexerTracks...

diff --git a/ITS/AliITSSecondaryVertices.C b/ITS/AliITSSecondaryVertices.C
new file mode 100644 (file)
index 0000000..b859208
--- /dev/null
+++ b/ITS/AliITSSecondaryVertices.C
@@ -0,0 +1,175 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+//-- --- standard headers------------- 
+#include <Riostream.h>
+//--------Root headers ---------------
+#include <TSystem.h>
+#include <TFile.h>
+#include <TStopwatch.h>
+#include <TObject.h>
+#include <TTree.h>
+#include <TChain.h>
+#include <TF1.h>
+#include <TLegend.h>
+#include <TLegendEntry.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TCut.h>
+#include <TParticle.h>
+#include <TCanvas.h>
+#include <TStyle.h>
+#include <TGraph.h>
+//----- AliRoot headers ---------------
+//#include "alles.h"
+#include "AliRun.h"
+#include "AliMagF.h"
+#include "AliKalmanTrack.h"
+#include "AliESDVertex.h"
+#include "AliITSVertexer.h"
+#include "AliITSVertexerTracks.h"
+#include <AliHeader.h>
+#include <AliGenEventHeader.h>
+#include "AliRun.h"
+#include "AliRunLoader.h"
+#include "AliITStrackV2.h"
+#include "AliITSSimpleVertex.h"
+#include "AliITSStrLine.h"
+#include "AliITSLoader.h"
+#include <AliESD.h>
+#include <AliStack.h>
+//-------------------------------------
+#endif  
+
+
+void AliITSSecondaryVertices(Int_t ntracks, Int_t *pos) {
+ 
+  // example of usige AliITSVertexerTracks to get a (secodnary) vertex from a set of tracks.
+
+  Double_t field = 0.4;
+ 
+  AliRunLoader *rl = AliRunLoader::Open("galice.root");
+  if (!rl) {
+    cerr<<"Can't start session !\n";
+  }
+
+  rl->LoadgAlice();
+  if (rl->GetAliRun()){
+    AliMagF * magf = gAlice->Field();
+    AliKalmanTrack:: SetFieldMap(magf);
+    AliKalmanTrack::SetUniformFieldTracking();
+  }else {
+    cerr<<"Can't get AliRun !\n";
+    
+  }
+
+  field=rl->GetAliRun()->Field()->SolenoidField()/10.;
+  rl->LoadHeader();
+ 
+  
+  TFile *ef=TFile::Open("AliESDs.root");
+  if (!ef || !ef->IsOpen()) {cerr<<"Can't AliESDs.root !n";}
+  AliESD* event = new AliESD;
+  TTree* tree = (TTree*) ef->Get("esdTree");
+  if (!tree) {cerr<<"no ESD tree foundn";};
+  tree->SetBranchAddress("ESD", &event);
+
+
+  Int_t e=0;
+  tree->GetEvent(e);
+  Int_t ntrk=event->GetNumberOfTracks();
+  cout<<"Number of ESD tracks : "<<ntrk<<endl; 
+
+  // Open input and output files
+  TFile *outFile = TFile::Open("AliITSVertices.root","recreate");
+    
+
+  TTree *trkTree = new TTree("TreeT_ITS","its tracks");
+  AliITStrackV2 *itstrack = 0;
+  trkTree->Branch("tracks","AliITStrackV2",&itstrack,ntrk,0);
+  Int_t pipe=0;
+  for(Int_t i=0; i<ntrk; i++) {
+    AliESDtrack *esdTrack = (AliESDtrack*)event->GetTrack(i);
+    UInt_t status=esdTrack->GetStatus();
+    UInt_t flags=AliESDtrack::kTPCin|AliESDtrack::kITSin;
+    if ((status&AliESDtrack::kTPCin)==0)continue;
+    if ((status&AliESDtrack::kITSin)==0)continue;    
+    if ((status&AliESDtrack::kITSrefit)==0) continue;
+    if ((status&flags)==status) pipe=1;
+
+    itstrack = new AliITStrackV2(*esdTrack);
+    if (pipe!=0) {
+      itstrack->PropagateTo(3.,0.0028,65.19);
+      itstrack->PropagateToVertex();  // This is not "exactly" the vertex 
+    }
+    Int_t nclus=itstrack->GetNumberOfClusters();
+
+    if(nclus<6)continue;
+    trkTree->Fill();
+    itstrack = 0;
+    // delete esdTrack; 
+  }
+  // Primary vertex
+  Double_t vtx[3];
+  event->GetVertex()->GetXYZ(vtx);
+  cout<<"Primary Vertex:"<<endl;
+  event->GetVertex()->PrintStatus();
+  cout<<"\n\n\n"<<endl;
+  
+
+  // Create vertexer
+  AliITSVertexerTracks *vertexer = new AliITSVertexerTracks();
+  vertexer->SetDCAcut(0);
+  vertexer->SetDebug(0);
+
+
+  Double_t vertF1[3],vertF2[3],vertF3[3],vertF4[3],vertF5[3];
+  Int_t ncombi1,ncombi2,ncombi3,ncombi4,ncombi5;
+  Double_t sig1,sig2,sig3,sig4,sig5; 
+  AliITSSimpleVertex *vert=0;
+
+  
+  // Calculate vertex from tracks in pos
+  // BEWARE: the method VertexForSelectedTracks returns the address of the data member fSimpVert 
+  // which is always the same for all calls to VertexForSelectedTracks.
+
+  cout<<"\nStraight Line Min Dist finder with weights"<<endl;
+  vert=(AliITSSimpleVertex*) vertexer->VertexForSelectedTracks(event,ntracks,pos,1);
+  vert->Print();
+  vert->GetXYZ(vertF1);
+  ncombi1=vert->GetNContributors();
+  sig1=vert->GetDispersion();
+
+
+  cout<<"Straight Line Min Dist finder"<<endl;
+  vert=(AliITSSimpleVertex*) vertexer->VertexForSelectedTracks(event,ntracks,pos,2);   
+  vert->Print();
+  vert->GetXYZ(vertF2);
+  ncombi2=vert->GetNContributors();
+  sig2=vert->GetDispersion();
+  
+  cout<<"Helix finder"<<endl;
+  vert=(AliITSSimpleVertex*) vertexer->VertexForSelectedTracks(event,ntracks,pos,3);   
+  vert->Print();
+  vert->GetXYZ(vertF3);
+  ncombi3=vert->GetNContributors();
+  sig3=vert->GetDispersion();
+ 
+  cout<<"Straight Line finder with weights"<<endl;
+  vert=(AliITSSimpleVertex*) vertexer->VertexForSelectedTracks(event,ntracks,pos,4);   
+  vert->Print();
+  vert->GetXYZ(vertF4);
+  ncombi4=vert->GetNContributors();
+  sig4=vert->GetDispersion();
+  
+  cout<<"Straight Line finder with weights"<<endl;
+  vert=(AliITSSimpleVertex*) vertexer->VertexForSelectedTracks(event,ntracks,pos,5);
+  vert->Print();
+  vert->GetXYZ(vertF5);
+  ncombi5=vert->GetNContributors();
+  sig5=vert->GetDispersion();
+  
+
+  delete vertexer;
+
+  outFile->Close();
+  return;
+}

diff --git a/ITS/AliITSSimpleVertex.cxx b/ITS/AliITSSimpleVertex.cxx
new file mode 100644 (file)
index 0000000..2e1d588
--- /dev/null
+++ b/ITS/AliITSSimpleVertex.cxx
@@ -0,0 +1,89 @@
+/**************************************************************************
+ * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+//-----------------------------------------------------------------
+//           Implementation of the Secondary Vertex class
+//           This class contains the Secondary Vertex
+//           of a set of tracks
+// Origin: F.Prino, Torino, prino@to.infn.it
+//-----------------------------------------------------------------
+
+#include "AliITSSimpleVertex.h"
+
+
+ClassImp(AliITSSimpleVertex)
+
+//--------------------------------------------------------------------------
+AliITSSimpleVertex::AliITSSimpleVertex() {
+//
+// Default Constructor, set everything to 0
+//
+  for(Int_t k=0;k<3;k++) fPosition[k]   = 0;
+  fSigma = 0;
+  fNContributors=0;
+}
+
+//--------------------------------------------------------------------------
+AliITSSimpleVertex::AliITSSimpleVertex(Double_t position[3],Double_t dispersion,
+               Int_t nContributors) {
+
+
+  //
+  // Constructor for vertex Z from pixels
+  //
+
+  for(Int_t k=0;k<3;k++) fPosition[k]   = position[k];
+  fSigma         = dispersion;
+  fNContributors = nContributors;
+
+}
+
+
+//--------------------------------------------------------------------------
+AliITSSimpleVertex::~AliITSSimpleVertex() {
+//  
+// Default Destructor
+//
+
+}
+//--------------------------------------------------------------------------
+void AliITSSimpleVertex::GetXYZ(Double_t position[3]) const {
+//
+// Return position of the vertex in global frame
+//
+  position[0] = fPosition[0];
+  position[1] = fPosition[1];
+  position[2] = fPosition[2];
+
+  return;
+}
+//--------------------------------------------------------------------------
+void AliITSSimpleVertex::Print(Option_t* /*option*/) const {
+//
+// Print out information on all data members
+//
+  printf("Secondary vertex position:\n");
+  printf("   x = %f\n",fPosition[0]);
+  printf("   y = %f\n",fPosition[1]);
+  printf("   z = %f\n",fPosition[2]);
+  printf(" Dispersion = %f\n",fSigma);
+  printf(" # tracks = %d\n",fNContributors);
+
+  return;
+}
+
+
+
+

diff --git a/ITS/AliITSSimpleVertex.h b/ITS/AliITSSimpleVertex.h
new file mode 100644 (file)
index 0000000..71c8c9a
--- /dev/null
+++ b/ITS/AliITSSimpleVertex.h
@@ -0,0 +1,67 @@
+#ifndef ALIITSSIMPLEVERTEX_H
+#define ALIITSSIMPLEVERTEX_H
+/* Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+
+//-------------------------------------------------------
+//                    Secondary Vertex Class
+//          as calculated by AliITSVertexerTracks
+//   Origin: F. Prino, Torino, prino@to.infn.it
+//-------------------------------------------------------
+
+/*****************************************************************************
+ *                                                                           *
+ * This class deals with secondary vertices.                                 *
+ * AliITSSimpleVertex objects are created by the class AliITSVertexerTracks  *
+ *                                                                           *
+ *****************************************************************************/
+
+//---- Root headers -----
+#include <TNamed.h>
+
+class AliITSSimpleVertex : public TNamed {
+ 
+ public:
+ 
+  AliITSSimpleVertex();
+  AliITSSimpleVertex(Double_t position[3],Double_t dispersion,
+               Int_t nContributors);
+  virtual ~AliITSSimpleVertex();
+
+
+  void   SetXYZ(Double_t pos[3]) {for(Int_t j=0; j<3; j++) fPosition[j]=pos[j];}
+  void   SetXv(Double_t xVert) {fPosition[0]=xVert; }
+  void   SetYv(Double_t yVert) {fPosition[1]=yVert; }
+  void   SetZv(Double_t zVert) {fPosition[2]=zVert; }
+  void   SetDispersion(Double_t disp) { fSigma=disp; }
+  void   SetNContributors(Int_t nContr) {fNContributors=nContr; }
+
+  void     GetXYZ(Double_t position[3]) const;
+  Double_t GetXv() const { return fPosition[0]; }
+  Double_t GetYv() const { return fPosition[1]; }
+  Double_t GetZv() const { return fPosition[2]; }
+  Double_t GetDispersion() const { return fSigma; }
+  Int_t    GetNContributors() const { return fNContributors; }
+
+  void     Print(Option_t* option = "") const;
+
+
+ protected:
+
+  Double_t fPosition[3];  // vertex position
+  Double_t fSigma;  // track dispersion around found vertex
+  Int_t    fNContributors;  // # of tracklets/tracks used for the estimate 
+
+
+  ClassDef(AliITSSimpleVertex,1)  // Class for Primary Vertex
+};
+
+#endif
+
+
+
+
+
+
+    

diff --git a/ITS/AliITSStrLine.cxx b/ITS/AliITSStrLine.cxx
index 3d912d7c20ade24de0b014567b98cd723bbf30e5..91abe642bba206a41aca36f943f7225d0b3c71cf 100644 (file)
--- a/ITS/AliITSStrLine.cxx
+++ b/ITS/AliITSStrLine.cxx
@@ -19,11 +19,7 @@
 //                                                               //
 ///////////////////////////////////////////////////////////////////
 #include <Riostream.h>
-#include <TMath.h>
-#include <TObject.h>
-#include <TObjArray.h>
 #include <TTree.h>
-#include <TClonesArray.h>
 #include "AliITSStrLine.h"
 
 ClassImp(AliITSStrLine)
@@ -76,6 +72,19 @@ void AliITSStrLine::PrintStatus() const {
   cout <<" Debug flag: "<<fDebug<<endl;
 }
 
+//________________________________________________________
+Int_t AliITSStrLine::IsParallelTo(AliITSStrLine *line) const {
+  // returns 1 if lines are parallel, 0 if not paralel
+  Double_t cd2[3];
+  line->GetCd(cd2);
+  Double_t vecpx=fCd[1]*cd2[2]-fCd[2]*cd2[1];
+  if(vecpx!=0) return 0;
+  Double_t vecpy=-fCd[0]*cd2[2]+fCd[2]*cd2[0];
+  if(vecpy!=0) return 0;
+  Double_t vecpz=fCd[0]*cd2[1]-fCd[1]*cd2[0];
+  if(vecpz!=0) return 0;
+  return 1;
+}
 //________________________________________________________
 Int_t AliITSStrLine::Crossrphi(AliITSStrLine *line){
   // Cross 2 lines in the X-Y plane
@@ -101,7 +110,7 @@ Int_t AliITSStrLine::Crossrphi(AliITSStrLine *line){
 }
 
 //________________________________________________________
-Int_t AliITSStrLine::Cross(AliITSStrLine *line, Double_t *point){
+Int_t AliITSStrLine::CrossPoints(AliITSStrLine *line, Double_t *point1, Double_t *point2){
   // Looks for the crossing point estimated starting from the
   // DCA segment
   Double_t p2[3];
@@ -125,15 +134,64 @@ Int_t AliITSStrLine::Cross(AliITSStrLine *line, Double_t *point){
   fTpar = (a11*k2-a21*k1) / deno;
   Double_t par2 = (k1*a22-k2*a12) / deno;
   line->SetPar(par2);
-  Double_t point1[3];
-  Double_t point2[3];
   GetCurrentPoint(point1);
   line->GetCurrentPoint(point2);
-  for(i=0;i<3;i++)point[i]=(point1[i]+point2[i])/2.;
   return 0;
 }
+//________________________________________________________________
+Int_t AliITSStrLine::Cross(AliITSStrLine *line, Double_t *point){
+
+  //Finds intersection between lines
+  Double_t point1[3];
+  Double_t point2[3];
+  Int_t retcod=CrossPoints(line,point1,point2);
+  if(retcod==0){
+    for(Int_t i=0;i<3;i++)point[i]=(point1[i]+point2[i])/2.;
+    return 0;
+  }else{
+    return -1;
+  }
+}
+
+//___________________________________________________________
+Double_t AliITSStrLine::GetDCA(AliITSStrLine *line){
+  //Returns the distance of closest approach between two lines
+  Double_t p2[3];
+  Double_t cd2[3];
+  line->GetP0(p2);
+  line->GetCd(cd2);
+  Int_t i;
+  Int_t ispar=IsParallelTo(line);
+  if(ispar){
+    Double_t dist1q=0,dist2=0,mod=0;
+    for(i=0;i<3;i++){
+      dist1q+=(fP0[i]-p2[i])*(fP0[i]-p2[i]);
+      dist2+=(fP0[i]-p2[i])*fCd[i];
+      mod+=fCd[i]*fCd[i];
+    }
+    if(mod!=0){
+      dist2/=mod;
+      return TMath::Sqrt(dist1q-dist2*dist2);
+    }else{return -1;}
+  }else{
+     Double_t perp[3];
+     perp[0]=fCd[1]*cd2[2]-fCd[2]*cd2[1];
+     perp[1]=-fCd[0]*cd2[2]+fCd[2]*cd2[0];
+     perp[2]=fCd[0]*cd2[1]-fCd[1]*cd2[0];
+     Double_t mod=0,dist=0;
+     for(i=0;i<3;i++){
+       mod+=perp[i]*perp[i];
+       dist+=(fP0[i]-p2[i])*perp[i];
+     }
+     mod=sqrt(mod);
+     if(mod!=0){
+       dist/=mod;
+       return TMath::Abs(dist);
+     }else{return -1;}
+  }
+}
 //________________________________________________________
-void AliITSStrLine::GetCurrentPoint(Double_t *point){
+void AliITSStrLine::GetCurrentPoint(Double_t *point) const {
   // Fills the array point with the current value on the line
   for(Int_t i=0;i<3;i++)point[i]=fP0[i]+fCd[i]*fTpar;
 }

diff --git a/ITS/AliITSStrLine.h b/ITS/AliITSStrLine.h
index 3654d6ebf90d4e095572d341ae0c0a425b3964d8..c0876f42f450a5f4ba5b37bb597b1855bbfdafcd 100644 (file)
--- a/ITS/AliITSStrLine.h
+++ b/ITS/AliITSStrLine.h
@@ -21,18 +21,22 @@ class AliITSStrLine : public TObject {
     void SetP0(Double_t *point) {for(Int_t i=0;i<3;i++)fP0[i]=point[i];}
     void SetCd(Double_t *cd) {for(Int_t i=0;i<3;i++)fCd[i]=cd[i];}
     void SetDebug(Int_t dbfl = 0){fDebug = dbfl; }  
-    void GetP0(Double_t *point) {for(Int_t i=0;i<3;i++)point[i]=fP0[i];}
-    void GetCd(Double_t *cd) {for(Int_t i=0;i<3;i++)cd[i]=fCd[i];}
-    void GetCurrentPoint(Double_t *point);
+    void GetP0(Double_t *point) const {for(Int_t i=0;i<3;i++)point[i]=fP0[i];}
+    void GetCd(Double_t *cd) const {for(Int_t i=0;i<3;i++)cd[i]=fCd[i];}
+    void GetCurrentPoint(Double_t *point) const;
+    Int_t IsParallelTo(AliITSStrLine *line) const;
     Int_t Crossrphi(AliITSStrLine *line);
-    Int_t Cross(AliITSStrLine *line,Double_t *point);
- private:
-    void SetPar(Double_t par){fTpar = par;}
+    Int_t CrossPoints(AliITSStrLine *line, Double_t *point1, Double_t *point2);
+    Int_t Cross(AliITSStrLine *line, Double_t *point);
+    Double_t GetDCA(AliITSStrLine *line);
  protected:
     Double_t fP0[3];           // given point
     Double_t fCd[3];           // direction cosines
     Double_t fTpar;            //! parameter 
     Int_t   fDebug;           //! debug flag - verbose printing if >0
+ private:
+    void SetPar(Double_t par){fTpar = par;}
+
   ClassDef(AliITSStrLine,1);
 };
 

diff --git a/ITS/AliITSVertexerTracks.cxx b/ITS/AliITSVertexerTracks.cxx
index 5b772cfbb375b246f83e47b43db139fbf79ebc30..7816e6b1f12ed86c61bd1a633679cb9479cbdff8 100644 (file)
--- a/ITS/AliITSVertexerTracks.cxx
+++ b/ITS/AliITSVertexerTracks.cxx
@@ -53,7 +53,7 @@ AliITSVertexerTracks::AliITSVertexerTracks():AliITSVertexer() {
   SetMinTracks();
   fTrksToSkip = 0;
   fNTrksToSkip = 0;
-  for(Int_t i=0; i<3; i++)fInitPos[i] = 0.;
+  fDCAcut=0;
 }
 //----------------------------------------------------------------------------
 AliITSVertexerTracks::AliITSVertexerTracks(TFile *inFile,TFile *outFile,
@@ -73,7 +73,7 @@ AliITSVertexerTracks::AliITSVertexerTracks(TFile *inFile,TFile *outFile,
   SetMinTracks();
   fTrksToSkip = 0;
   fNTrksToSkip = 0;
-  for(Int_t i=0; i<3; i++) fInitPos[i] = 0.;
+  fDCAcut=0;
   SetDebug();
 }
 //----------------------------------------------------------------------------
@@ -90,7 +90,7 @@ AliITSVertexerTracks::AliITSVertexerTracks(const AliMagF *map, TString fn,
   SetMinTracks();
   fTrksToSkip = 0;
   fNTrksToSkip = 0;
-  for(Int_t i=0; i<3; i++) fInitPos[i] = 0.;
+  fDCAcut=0;
 }
 //______________________________________________________________________
 AliITSVertexerTracks::AliITSVertexerTracks(const AliITSVertexerTracks &vtxr) : AliITSVertexer(vtxr) {
@@ -158,7 +158,7 @@ void AliITSVertexerTracks::FindVertices() {
   for(Int_t ev=fFirstEvent; ev<=fLastEvent; ev++) {
     if(ev % 100 == 0 || fDebug) printf("--- Processing event %d of %d ---\n",ev,fLastEvent);
 
-    FindVertexForCurrentEvent(ev);
+    FindPrimaryVertexForCurrentEvent(ev);
 
     if(!fCurrentVertex) {
       printf("AliITSVertexerTracks::FindVertixes(): no tracks tree for event %d\n",ev);
@@ -208,14 +208,14 @@ void AliITSVertexerTracks::FindVerticesESD() {
     if(ev % 100 == 0 || fDebug) 
       printf("--- Processing event %d of %d ---\n",ev,fLastEvent-fFirstEvent);
 
-    FindVertexForCurrentEvent(esdEvent);
+    FindPrimaryVertexForCurrentEvent(esdEvent);
 
     if(!fCurrentVertex) {
       printf("AliITSVertexerTracks::FindVertixesESD():\n no vertex for event %d\n",ev);
       continue;
     }
 
-    cout<<"VERTICE TROVATO\n";
+
     if(fDebug) fCurrentVertex->PrintStatus();
 
     // write the ESD to file
@@ -233,14 +233,13 @@ void AliITSVertexerTracks::FindVerticesESD() {
 }
 //----------------------------------------------------------------------------
 Int_t AliITSVertexerTracks::PrepareTracks(TTree &trkTree) {
-//
-// Propagate tracks to initial vertex position and store them in a TObjArray
-//
-  Double_t maxd0rphi = 3.;
-  Double_t alpha,xlStart,d0rphi;
+  //
+  // Propagate tracks to initial vertex position and store them in a TObjArray
+  //
+  Double_t maxd0rphi = 3.;  
   Int_t    nTrks    = 0;
   Bool_t   skipThis;
-
+  Double_t d0rphi;
   Int_t    nEntries = (Int_t)trkTree.GetEntries();
 
   if(!fTrkArray.IsEmpty()) fTrkArray.Clear();
@@ -250,7 +249,7 @@ Int_t AliITSVertexerTracks::PrepareTracks(TTree &trkTree) {
     printf(" PrepareTracks()\n");
     trkTree.Print();
   }
-
+  cout<<" entr tree its tracks = "<<nEntries<<endl;
   for(Int_t i=0; i<nEntries; i++) {
     // check tracks to skip
     skipThis = kFALSE;
@@ -265,41 +264,89 @@ Int_t AliITSVertexerTracks::PrepareTracks(TTree &trkTree) {
     AliITStrackV2 *itstrack = new AliITStrackV2; 
     trkTree.SetBranchAddress("tracks",&itstrack);
     trkTree.GetEvent(i);
+    d0rphi=Prepare(itstrack);
+    if(d0rphi> maxd0rphi) { if(fDebug)cout<<"    !!!! d0rphi "<<d0rphi<<endl;continue; }
+    fTrkArray.AddLast(itstrack);
+    
+    nTrks++; 
+    if(fDebug)cout<<" :-) nTrks, d0rphi "<<nTrks<<"  "<<d0rphi<<endl;
+   
+  }
+  if(fTrksToSkip) delete [] fTrksToSkip;
 
+  return nTrks;
+} 
+//----------------------------------------------------------------------------
+Int_t AliITSVertexerTracks::PrepareTracks(AliESD* esdEvent,Int_t nofCand, Int_t *trkPos) {
+  //
+  // Propagate tracks to initial vertex position and store them in a TObjArray
+  //
+  Int_t    nTrks    = 0;
+  Double_t maxd0rphi = 3.; 
+  Double_t d0rphi; 
 
-    // propagate track to vtxSeed
-    alpha  = itstrack->GetAlpha();
-    xlStart = fNominalPos[0]*TMath::Cos(alpha)+fNominalPos[1]*TMath::Sin(alpha);
-    itstrack->PropagateTo(3.,0.0023,65.19); // to beam pipe (0.8 mm of Be) 
-    itstrack->PropagateTo(xlStart,0.,0.);   // to vtxSeed
+  if(!fTrkArray.IsEmpty()) fTrkArray.Clear();
+  fTrkArray.Expand(100);
 
-    // select tracks with d0rphi < maxd0rphi
-    d0rphi = TMath::Abs(itstrack->GetD(fNominalPos[0],fNominalPos[1]));
-    if(d0rphi > maxd0rphi) { delete itstrack; continue; }
-   
+  if(fDebug) {
+    printf(" PrepareTracks()\n");
+  }
+  AliITStrackV2* itstrack;
+  
+  for(Int_t i=0; i<nofCand;i++){
+    AliESDtrack *esdTrack = (AliESDtrack*)esdEvent->GetTrack(trkPos[i]);
+    UInt_t status=esdTrack->GetStatus();
+    if ((status&AliESDtrack::kTPCin)==0)continue;
+    if ((status&AliESDtrack::kITSin)==0)continue;
+    if ((status&AliESDtrack::kITSrefit)==0) continue;
+
+    itstrack = new AliITStrackV2(*esdTrack);
+    d0rphi=Prepare(itstrack);
+    if(d0rphi> maxd0rphi) { if(fDebug)cout<<"    !!!! d0rphi "<<d0rphi<<endl;continue; }
+    Int_t nclus=itstrack->GetNumberOfClusters();
+
+    if(nclus<6){delete itstrack;continue;}
     fTrkArray.AddLast(itstrack);
-
+    
     nTrks++; 
+    if(fDebug)cout<<" :-) nTrks, d0rphi "<<nTrks<<"  "<<d0rphi<<endl;
+    //delete itstrack;
   }
+ 
 
   if(fTrksToSkip) delete [] fTrksToSkip;
-
   return nTrks;
-} 
+}
+//----------------------------------------------------------------------------
+Double_t AliITSVertexerTracks::Prepare(AliITStrackV2* itstrack){
+ 
+  //
+  Double_t alpha,xlStart,d0rphi; 
+  // propagate track to vtxSeed
+  alpha  = itstrack->GetAlpha();
+  xlStart = fNominalPos[0]*TMath::Cos(alpha)+fNominalPos[1]*TMath::Sin(alpha);
+  if(itstrack->GetX()>3.)itstrack->PropagateTo(3.,0.0023,65.19); // to beam pipe (0.8 mm of Be) 
+  itstrack->PropagateTo(xlStart,0.,0.); 
+  // select tracks with d0rphi < maxd0rphi
+  d0rphi = TMath::Abs(itstrack->GetD(fNominalPos[0],fNominalPos[1]));
+  return d0rphi;
+          
+}
 //----------------------------------------------------------------------------
 void AliITSVertexerTracks::PrintStatus() const {
 //
 // Print status
 //
   printf(" Initial position (%f,%f)\n",fNominalPos[0],fNominalPos[1]);
-  printf(" Vertex position after vertex finder (%f, %f, %f)\n",fInitPos[0],fInitPos[1],fInitPos[2]);
+  printf(" Vertex position after vertex finder:\n");
+  fSimpVert.Print();
   printf(" Number of tracks in array: %d\n",(Int_t)fTrkArray.GetEntriesFast());
   printf(" Minimum # tracks required in fit: %d\n",fMinTracks);
 
   return;
 }
 //----------------------------------------------------------------------------
-AliESDVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(Int_t evnumb) {
+AliESDVertex* AliITSVertexerTracks::FindPrimaryVertexForCurrentEvent(Int_t evnumb) {
 //
 // Vertex for current event
 //
@@ -329,7 +376,7 @@ AliESDVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(Int_t evnumb) {
   return fCurrentVertex;
 }
 //----------------------------------------------------------------------------
-AliESDVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(AliESD *esdEvent)
+AliESDVertex* AliITSVertexerTracks::FindPrimaryVertexForCurrentEvent(AliESD *esdEvent)
 {
 //
 // Vertex for current ESD event
@@ -351,7 +398,7 @@ AliESDVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(AliESD *esdEvent)
       itstrack = new AliITStrackV2(*esdTrack);
     }
     catch (const Char_t *msg) {
-        Warning("FindVertexForCurrentEvent",msg);
+        Warning("FindPrimaryVertexForCurrentEvent",msg);
         delete esdTrack;
         continue;
     }
@@ -371,7 +418,6 @@ AliESDVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(AliESD *esdEvent)
   // Set initial vertex position from ESD
   esdEvent->GetVertex()->GetXYZ(vtx);
   SetVtxStart(vtx[0],vtx[1]);
-
   // VERTEX FINDER
   VertexFinder();
 
@@ -393,6 +439,37 @@ AliESDVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(AliESD *esdEvent)
   cout<<"Vertex: "<<vtx[0]<<", "<<vtx[1]<<", "<<vtx[2]<<endl;
   return fCurrentVertex;
 }
+//----------------------------------------------------------------------------
+AliITSSimpleVertex* AliITSVertexerTracks::VertexForSelectedTracks(AliESD *esdEvent,Int_t nofCand, Int_t *trkPos,  Int_t opt){
+
+  //
+  // Computes the vertex for selected tracks 
+  // trkPos=vector with track positions in ESD
+  // values of opt -> see AliITSVertexerTracks.h
+  //
+  Double_t vtx[3]={0,0,0};
+
+  Int_t nTrks = PrepareTracks(esdEvent,nofCand, trkPos);
+  //delete trkTree;//  :-)) 
+  if(fDebug) printf(" tracks prepared: %d\n",nTrks);
+  if(nTrks < fMinTracks) {
+    fSimpVert.SetXYZ(vtx);
+    fSimpVert.SetDispersion(999);
+    fSimpVert.SetNContributors(-5);
+    return &fSimpVert;
+  }
+ 
+  // Set initial vertex position from ESD
+  esdEvent->GetVertex()->GetXYZ(vtx);
+  SetVtxStart(vtx[0],vtx[1]);
+  if(opt==1)  StrLinVertexFinderMinDist(1);
+  if(opt==2)  StrLinVertexFinderMinDist(0);
+  if(opt==3)  HelixVertexFinder();
+  if(opt==4)  VertexFinder(1);
+  if(opt==5)  VertexFinder(0);
+  return &fSimpVert;
+}
+
 //---------------------------------------------------------------------------
 void AliITSVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) {
 //
@@ -413,34 +490,33 @@ void AliITSVertexerTracks::TooFewTracks() {
   return;
 }
 //---------------------------------------------------------------------------
-void AliITSVertexerTracks::VertexFinder() {
+void AliITSVertexerTracks::VertexFinder(Int_t OptUseWeights) {
 
   // Get estimate of vertex position in (x,y) from tracks DCA
-  // Then this estimate is stored to the data member fInitPos   
+  // Then this estimate is stored to the data member fSimpVert  
   // (previous values are overwritten)
 
-
  
-  /*
-******* TEMPORARY!!! FOR TEST ONLY!!! **********************************
-
-fInitPos[0] = fNominalPos[0]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing
-fInitPos[1] = fNominalPos[1]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing
-  */
-
-  fInitPos[2] = 0.;
-  for(Int_t i=0;i<2;i++)fInitPos[i]=fNominalPos[i];
-
+  Double_t initPos[3];
+  initPos[2] = 0.;
+  for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
   Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
-
   Double_t aver[3]={0.,0.,0.};
+  Double_t aversq[3]={0.,0.,0.};
+  Double_t sigmasq[3]={0.,0.,0.};
+  Double_t sigma=0;
   Int_t ncombi = 0;
   AliITStrackV2 *track1;
   AliITStrackV2 *track2;
-  AliITSStrLine line1;
-  AliITSStrLine line2;
+  Double_t alpha,mindist;
+
   for(Int_t i=0; i<nacc; i++){
     track1 = (AliITStrackV2*)fTrkArray.At(i);
+    alpha=track1->GetAlpha();
+    mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
+    AliITSStrLine *line1 = new AliITSStrLine();
+    track1->ApproximateHelixWithLine(mindist,line1);
+   
     if(fDebug>5){
       Double_t xv,par[5];
       track1->GetExternalParameters(xv,par);
@@ -448,60 +524,335 @@ fInitPos[1] = fNominalPos[1]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian sm
       for(Int_t ii=0;ii<5;ii++)cout<<par[ii]<<" ";
       cout<<endl;
     }
-    Double_t mom1[3];
-    Double_t alpha = track1->GetAlpha();
-    Double_t azim = TMath::ASin(track1->GetSnp())+alpha;
-    Double_t theta = TMath::Pi()/2. - TMath::ATan(track1->GetTgl());
-    mom1[0] = TMath::Sin(theta)*TMath::Cos(azim);
-    mom1[1] = TMath::Sin(theta)*TMath::Sin(azim);
-    mom1[2] = TMath::Cos(theta);
-
-    Double_t pos1[3];
-    Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
-    track1->GetGlobalXYZat(mindist,pos1[0],pos1[1],pos1[2]);
-    line1.SetP0(pos1);
-    line1.SetCd(mom1);
+
     for(Int_t j=i+1; j<nacc; j++){
       track2 = (AliITStrackV2*)fTrkArray.At(j);
-      Double_t mom2[3];
-      alpha = track2->GetAlpha();
-      azim = TMath::ASin(track2->GetSnp())+alpha;
-      theta = TMath::Pi()/2. - TMath::ATan(track2->GetTgl());
-      mom2[0] = TMath::Sin(theta)*TMath::Cos(azim);
-      mom2[1] = TMath::Sin(theta)*TMath::Sin(azim);
-      mom2[2] = TMath::Cos(theta);
-      Double_t pos2[3];
+      alpha=track2->GetAlpha();
       mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
-      track2->GetGlobalXYZat(mindist,pos2[0],pos2[1],pos2[2]);
-      line2.SetP0(pos2);
-      line2.SetCd(mom2);
-      Double_t crosspoint[3];
-      Int_t retcode = line2.Cross(&line1,crosspoint);
-      if(retcode<0){
-       if(fDebug>10)cout<<" i= "<<i<<",   j= "<<j<<endl;
-       if(fDebug>10)cout<<"bad intersection\n";
-       line1.PrintStatus();
-       line2.PrintStatus();
+      AliITSStrLine *line2 = new AliITSStrLine();
+      track2->ApproximateHelixWithLine(mindist,line2);
+      Double_t distCA=line2->GetDCA(line1);
+      if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
+       Double_t pnt1[3],pnt2[3],crosspoint[3];
+
+       if(OptUseWeights<=0){
+         Int_t retcode = line2->Cross(line1,crosspoint);
+         if(retcode<0){
+           if(fDebug>10)cout<<" i= "<<i<<",   j= "<<j<<endl;
+           if(fDebug>10)cout<<"bad intersection\n";
+           line1->PrintStatus();
+           line2->PrintStatus();
+         }
+         else {
+           ncombi++;
+           for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
+           for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
+           if(fDebug>10)cout<<" i= "<<i<<",   j= "<<j<<endl;
+           if(fDebug>10)cout<<"\n Cross point: ";
+           if(fDebug>10)cout<<crosspoint[0]<<" "<<crosspoint[1]<<" "<<crosspoint[2]<<endl;
+         }
+       }
+       if(OptUseWeights>0){
+         Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2);
+         if(retcode>=0){
+           Double_t alpha, cs, sn;
+           alpha=track1->GetAlpha();
+           cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);   
+           Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
+           alpha=track2->GetAlpha();
+           cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
+           Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
+           Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
+           Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
+           Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
+           Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
+           crosspoint[0]=wx1*pnt1[0] + wx2*pnt2[0]; 
+           crosspoint[1]=wy1*pnt1[1] + wy2*pnt2[1]; 
+           crosspoint[2]=wz1*pnt1[2] + wz2*pnt2[2];
+         
+           ncombi++;
+           for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
+           for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
+         }
+       }
       }
-      else {
+      delete line2;
+    }
+    delete line1;
+  }
+  if(ncombi>0){
+    for(Int_t jj=0;jj<3;jj++){
+      initPos[jj] = aver[jj]/ncombi;
+      aversq[jj]/=ncombi;
+      sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
+      sigma+=sigmasq[jj];
+    }
+    sigma=TMath::Sqrt(TMath::Abs(sigma));
+  }
+  else {
+    Warning("VertexFinder","Finder did not succed");
+    sigma=999;
+  }
+  fSimpVert.SetXYZ(initPos);
+  fSimpVert.SetDispersion(sigma);
+  fSimpVert.SetNContributors(ncombi);
+}
+//---------------------------------------------------------------------------
+void AliITSVertexerTracks::HelixVertexFinder() {
+
+  // Get estimate of vertex position in (x,y) from tracks DCA
+  // Then this estimate is stored to the data member fSimpVert  
+  // (previous values are overwritten)
+
+
+  Double_t initPos[3];
+  initPos[2] = 0.;
+  for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
+
+  Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
+
+  Double_t aver[3]={0.,0.,0.};
+  Double_t averquad[3]={0.,0.,0.};
+  Double_t sigmaquad[3]={0.,0.,0.};
+  Double_t sigma=0;
+  Int_t ncombi = 0;
+  AliITStrackV2 *track1;
+  AliITStrackV2 *track2;
+  Double_t distCA;
+  Double_t x, par[5];
+  Double_t alpha, cs, sn;
+  Double_t crosspoint[3];
+  for(Int_t i=0; i<nacc; i++){
+    track1 = (AliITStrackV2*)fTrkArray.At(i);
+    
+
+    for(Int_t j=i+1; j<nacc; j++){
+      track2 = (AliITStrackV2*)fTrkArray.At(j);
+      
+
+      distCA=track2->PropagateToDCA(track1);
+
+      if(fDCAcut<=0 ||(fDCAcut>0&&distCA<fDCAcut)){
+       track1->GetExternalParameters(x,par);
+       alpha=track1->GetAlpha();
+       cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
+       Double_t x1=x*cs - par[0]*sn;
+       Double_t y1=x*sn + par[0]*cs;
+       Double_t z1=par[1];
+       Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2(); 
+
+       track2->GetExternalParameters(x,par);
+       alpha=track2->GetAlpha();
+       cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
+       Double_t x2=x*cs - par[0]*sn;
+       Double_t y2=x*sn + par[0]*cs;
+       Double_t z2=par[1];
+       Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2(); 
+       //      printf("Track %d pos=(%f,%f,%f) - dca=%f\n",i,x1,y1,z1,distCA);
+       //printf("Track %d pos=(%f,%f,%f)\n",j,x2,y2,z2);
+
+       Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
+       Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
+       Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
+       Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
+       crosspoint[0]=wx1*x1 + wx2*x2; 
+       crosspoint[1]=wy1*y1 + wy2*y2; 
+       crosspoint[2]=wz1*z1 + wz2*z2;
+
        ncombi++;
        for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
-       if(fDebug>10)cout<<" i= "<<i<<",   j= "<<j<<endl;
-       if(fDebug>10)cout<<"\n Cross point: ";
-       if(fDebug>10)cout<<crosspoint[0]<<" "<<crosspoint[1]<<" "<<crosspoint[2]<<endl;
+       for(Int_t jj=0;jj<3;jj++)averquad[jj]+=(crosspoint[jj]*crosspoint[jj]);
       }
     }
+      
   }
   if(ncombi>0){
-    for(Int_t jj=0;jj<3;jj++)fInitPos[jj] = aver[jj]/ncombi;
+    for(Int_t jj=0;jj<3;jj++){
+      initPos[jj] = aver[jj]/ncombi;
+      averquad[jj]/=ncombi;
+      sigmaquad[jj]=averquad[jj]-initPos[jj]*initPos[jj];
+      sigma+=sigmaquad[jj];
+    }
+    sigma=TMath::Sqrt(TMath::Abs(sigma));
   }
   else {
     Warning("VertexFinder","Finder did not succed");
+    sigma=999;
   }
+  fSimpVert.SetXYZ(initPos);
+  fSimpVert.SetDispersion(sigma);
+  fSimpVert.SetNContributors(ncombi);
+}
+//---------------------------------------------------------------------------
+void AliITSVertexerTracks::StrLinVertexFinderMinDist(Int_t OptUseWeights){
 
+  // Calculate the point at minimum distance to prepared tracks 
+  // Then this estimate is stored to the data member fSimpVert  
+  // (previous values are overwritten)
+  
+  Double_t initPos[3];
+  initPos[2] = 0.;
+  Double_t sigma=0;
+  for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
+  const Int_t knacc = (Int_t)fTrkArray.GetEntriesFast();
 
-  //************************************************************************
-  return;
+  AliITStrackV2 *track1;
+  Double_t (*vectP0)[3]=new Double_t [knacc][3];
+  Double_t (*vectP1)[3]=new Double_t [knacc][3];
+  
+  Double_t sum[3][3];
+  Double_t dsum[3]={0,0,0};
+  for(Int_t i=0;i<3;i++)
+    for(Int_t j=0;j<3;j++)sum[i][j]=0;
+  for(Int_t i=0; i<knacc; i++){
+    track1 = (AliITStrackV2*)fTrkArray.At(i);
+    Double_t alpha=track1->GetAlpha();
+    Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
+    AliITSStrLine *line1 = new AliITSStrLine();
+    track1->ApproximateHelixWithLine(mindist,line1);
+
+    Double_t p0[3],cd[3];
+    line1->GetP0(p0);
+    line1->GetCd(cd);
+    Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
+    vectP0[i][0]=p0[0];
+    vectP0[i][1]=p0[1];
+    vectP0[i][2]=p0[2];
+    vectP1[i][0]=p1[0];
+    vectP1[i][1]=p1[1];
+    vectP1[i][2]=p1[2];
+    
+    Double_t matr[3][3];
+    Double_t dknow[3];
+    if(OptUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
+    if(OptUseWeights==1){
+      Double_t sigmasq[3];
+      sigmasq[0]=track1->GetSigmaY2();
+      sigmasq[1]=track1->GetSigmaY2();
+      sigmasq[2]=track1->GetSigmaZ2();
+      GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
+    }
+
+    for(Int_t iii=0;iii<3;iii++){
+      dsum[iii]+=dknow[iii]; 
+      for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
+    }
+    delete line1;
+  }
+  
+  Double_t vett[3][3];
+  Double_t det=GetDeterminant3X3(sum);
+  
+   if(det!=0){
+     for(Int_t zz=0;zz<3;zz++){
+       for(Int_t ww=0;ww<3;ww++){
+        for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
+       }
+       for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
+       initPos[zz]=GetDeterminant3X3(vett)/det;
+     }
+
+
+     for(Int_t i=0; i<knacc; i++){
+       Double_t p0[3]={0,0,0},p1[3]={0,0,0};
+       for(Int_t ii=0;ii<3;ii++){
+        p0[ii]=vectP0[i][ii];
+        p1[ii]=vectP1[i][ii];
+       }
+       sigma+=GetStrLinMinDist(p0,p1,initPos);
+     }
+
+     sigma=TMath::Sqrt(sigma);
+   }else{
+    Warning("VertexFinder","Finder did not succed");
+    sigma=999;
+  }
+  delete vectP0;
+  delete vectP1;
+  fSimpVert.SetXYZ(initPos);
+  fSimpVert.SetDispersion(sigma);
+  fSimpVert.SetNContributors(knacc);
+}
+//_______________________________________________________________________
+Double_t AliITSVertexerTracks::GetDeterminant3X3(Double_t matr[][3]){
+  //
+  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];
+ return det;
+}
+//____________________________________________________________________________
+void AliITSVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t m[][3],Double_t *d){
+
+  //
+  Double_t x12=p0[0]-p1[0];
+  Double_t y12=p0[1]-p1[1];
+  Double_t z12=p0[2]-p1[2];
+  Double_t kk=x12*x12+y12*y12+z12*z12;
+  m[0][0]=2-2/kk*x12*x12;
+  m[0][1]=-2/kk*x12*y12;
+  m[0][2]=-2/kk*x12*z12;
+  m[1][0]=-2/kk*x12*y12;
+  m[1][1]=2-2/kk*y12*y12;
+  m[1][2]=-2/kk*y12*z12;
+  m[2][0]=-2/kk*x12*z12;
+  m[2][1]=-2*y12*z12;
+  m[2][2]=2-2/kk*z12*z12;
+  d[0]=2*p0[0]-2/kk*p0[0]*x12*x12-2/kk*p0[2]*x12*z12-2/kk*p0[1]*x12*y12;
+  d[1]=2*p0[1]-2/kk*p0[1]*y12*y12-2/kk*p0[0]*x12*y12-2/kk*p0[2]*z12*y12;
+  d[2]=2*p0[2]-2/kk*p0[2]*z12*z12-2/kk*p0[0]*x12*z12-2/kk*p0[1]*z12*y12;
+
+}
+//____________________________________________________________________________
+void AliITSVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t *sigmasq,Double_t m[][3],Double_t *d){
+  //
+  Double_t x12=p1[0]-p0[0];
+  Double_t y12=p1[1]-p0[1];
+  Double_t z12=p1[2]-p0[2];
+
+  Double_t den= x12*x12*sigmasq[1]*sigmasq[2]+y12*y12*sigmasq[0]*sigmasq[2]+z12*z12*sigmasq[0]*sigmasq[1];
+
+  Double_t kk= 2*(x12*x12/sigmasq[0]+y12*y12/sigmasq[1]+z12*z12/sigmasq[2]);
+
+  Double_t cc[3];
+  cc[0]=-x12/sigmasq[0];
+  cc[1]=-y12/sigmasq[1];
+  cc[2]=-z12/sigmasq[2];
+
+  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;
+
+  Double_t ss= -p0[0]*cc[0]-p0[1]*cc[1]-p0[2]*cc[2];
+
+  Double_t aa[3];
+  aa[0]=x12*sigmasq[1]*sigmasq[2]/den;
+  aa[1]=y12*sigmasq[0]*sigmasq[2]/den;
+  aa[2]=z12*sigmasq[0]*sigmasq[1]/den;
+
+  m[0][0]=aa[0]*(aa[0]*kk+2*cc[0])+2*cc[0]*aa[0]+2/sigmasq[0];
+  m[0][1]=aa[1]*(aa[0]*kk+2*cc[0])+2*cc[1]*aa[0];
+  m[0][2]=aa[2]*(aa[0]*kk+2*cc[0])+2*cc[2]*aa[0];
+
+  m[1][0]=aa[0]*(aa[1]*kk+2*cc[1])+2*cc[0]*aa[1];
+  m[1][1]=aa[1]*(aa[1]*kk+2*cc[1])+2*cc[1]*aa[1]+2/sigmasq[1];
+  m[1][2]=aa[2]*(aa[1]*kk+2*cc[1])+2*cc[2]*aa[1];
+
+  m[2][0]=aa[0]*(aa[2]*kk+2*cc[2])+2*cc[0]*aa[2];
+  m[2][1]=aa[1]*(aa[2]*kk+2*cc[2])+2*cc[1]*aa[2];
+  m[2][2]=aa[2]*(aa[2]*kk+2*cc[2])+2*cc[2]*aa[2]+2/sigmasq[2];
+
+  d[0]=-ww*(aa[0]*kk+2*cc[0])-2*ss*aa[0]+2*p0[0]/sigmasq[0];
+  d[1]=-ww*(aa[1]*kk+2*cc[1])-2*ss*aa[1]+2*p0[1]/sigmasq[1];
+  d[2]=-ww*(aa[2]*kk+2*cc[2])-2*ss*aa[2]+2*p0[2]/sigmasq[2];
+
+}
+//_____________________________________________________________________________
+Double_t AliITSVertexerTracks::GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t *x0){
+  //
+  Double_t x12=p0[0]-p1[0];
+  Double_t y12=p0[1]-p1[1];
+  Double_t z12=p0[2]-p1[2];
+  Double_t x10=p0[0]-x0[0];
+  Double_t y10=p0[1]-x0[1];
+  Double_t z10=p0[2]-x0[2];
+  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);
 
 }
 //---------------------------------------------------------------------------
@@ -520,8 +871,10 @@ void AliITSVertexerTracks::VertexFitter() {
   Int_t i,j,k,step=0;
   TMatrixD rv(3,1);
   TMatrixD vV(3,3);
-  rv(0,0) = fInitPos[0];
-  rv(1,0) = fInitPos[1];
+  Double_t initPos[3];
+  fSimpVert.GetXYZ(initPos);
+  rv(0,0) = initPos[0];
+  rv(1,0) = initPos[1];
   rv(2,0) = 0.;
   Double_t xlStart,alpha;
   Double_t rotAngle;
@@ -558,7 +911,7 @@ void AliITSVertexerTracks::VertexFitter() {
       // get track from track array
       t = (AliITStrackV2*)fTrkArray.At(k);
       alpha = t->GetAlpha();
-      xlStart = fInitPos[0]*TMath::Cos(alpha)+fInitPos[1]*TMath::Sin(alpha);
+      xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
       t->PropagateTo(xlStart,0.,0.);   // to vtxSeed
       rotAngle = alpha;
       if(alpha<0.) rotAngle += 2.*TMath::Pi();

diff --git a/ITS/AliITSVertexerTracks.h b/ITS/AliITSVertexerTracks.h
index f51939efe53d8c823d2d1e53bb3f422dad6f8115..3f8e5b61237f6e5a6aabad90b69b9ddcb7a28ca5 100644 (file)
--- a/ITS/AliITSVertexerTracks.h
+++ b/ITS/AliITSVertexerTracks.h
@@ -27,12 +27,15 @@
 
 #include "AliKalmanTrack.h"
 #include "AliITSVertexer.h"
+#include "AliITSSimpleVertex.h"
 
 #include <TObjArray.h>
 
 class TTree; 
 class AliESDVertex; 
+class AliITSSimpleVertex; 
 class AliESD;
+class AliITStrackV2;
 
 class AliITSVertexerTracks : public AliITSVertexer {
   
@@ -51,19 +54,47 @@ class AliITSVertexerTracks : public AliITSVertexer {
   // return vertex from the set of tracks in the tree
   AliESDVertex* VertexOnTheFly(TTree &trkTree);
   // computes the vertex for the current event
-  virtual AliESDVertex* FindVertexForCurrentEvent(Int_t evnumb);
+  virtual AliESDVertex* FindPrimaryVertexForCurrentEvent(Int_t evnumb);
+  virtual AliESDVertex* FindVertexForCurrentEvent(Int_t evnumb){
+    Warning(" FindVertexForCurrentEvent","Deprecated method use FindPrimaryVertexForCurrentEvent instead");
+    return FindPrimaryVertexForCurrentEvent(evnumb);
+  }
   // computes the vertex for the current event using the ESD
-  AliESDVertex*         FindVertexForCurrentEvent(AliESD *esdEvent);
+  AliESDVertex*         FindPrimaryVertexForCurrentEvent(AliESD *esdEvent);
+  AliESDVertex*         FindVertexForCurrentEvent(AliESD *esdEvent){
+    Warning(" FindVertexForCurrentEvent","Deprecated method use FindPrimaryVertexForCurrentEvent instead");
+    return FindPrimaryVertexForCurrentEvent(esdEvent);
+  }
   // computes the vertex for each event and stores it on file
   virtual void  FindVertices();
   // computes the vertex for each event and stores it in the ESD
   void FindVerticesESD();
+
+  // computes the vertex for selected tracks 
+  // (TrkPos=vector with track positions in ESD)
+  AliITSSimpleVertex* VertexForSelectedTracks(AliESD *esdEvent,Int_t nofCand, Int_t *trkPos, Int_t opt=1); 
+  // opt=1 (default) finds minimum-distance point among all the selected tracks
+  //       approximated as straight lines 
+  //       and uses errors on track parameters as weights
+  // opt=2 finds minimum-distance point among all the selected tracks
+  //       approximated as straight lines 
+  // opt=3 finds the average point among DCA points of all pairs of tracks
+  //       treated as helices
+  // opt=4 finds the average point among DCA points of all pairs of tracks
+  //       approximated as straight lines 
+  //       and uses errors on track parameters as weights
+  // opt=5 finds the average point among DCA points of all pairs of tracks
+  //       approximated as straight lines 
+
+
   virtual void  PrintStatus() const;
   void  SetFieldMap(const AliMagF *map)const{AliKalmanTrack::SetFieldMap(map);}
   void  SetMinTracks(Int_t n=2) { fMinTracks = n; return; }
   void  SetSkipTracks(Int_t n,Int_t *skipped); 
   void  SetVtxStart(Double_t x=0,Double_t y=0) 
     { fNominalPos[0]=x; fNominalPos[1]=y; return; }
+  void  SetDCAcut(Double_t maxdca)
+    { fDCAcut=maxdca; return;}
   
  private:
     // copy constructor (NO copy allowed: the constructor is protected
@@ -73,9 +104,10 @@ class AliITSVertexerTracks : public AliITSVertexer {
     AliITSVertexerTracks& operator=(const AliITSVertexerTracks& /* vtxr */);
   TFile    *fInFile;          // input file (with tracks)
   TFile    *fOutFile;         // output file for vertices
-  Double_t  fInitPos[3];      // vertex position after vertex finder
+  AliITSSimpleVertex  fSimpVert; // vertex after vertex finder
   Double_t  fNominalPos[2];   // initial knowledge on vertex position
   Int_t     fMinTracks;       // minimum number of tracks
+  Double_t  fDCAcut;          // maximum DCA between 2 tracks used for vertex
   Double_t  fMaxChi2PerTrack; // maximum contribition to the chi2 
   TObjArray fTrkArray;        // array with tracks to be processed
   Int_t     *fTrksToSkip;     // tracks to be skipped for find and fit 
@@ -84,8 +116,17 @@ class AliITSVertexerTracks : public AliITSVertexer {
   Bool_t   CheckField() const; 
   void     ComputeMaxChi2PerTrack(Int_t nTracks);
   Int_t    PrepareTracks(TTree &trkTree);
+  Int_t    PrepareTracks(AliESD* esdEvent,Int_t NofCand, Int_t *TrkPos);
+  Double_t Prepare(AliITStrackV2* itstrack);
   void     TooFewTracks();
-  void     VertexFinder();
+  void     VertexFinder(Int_t OptUseWeights=0);
+  void     HelixVertexFinder();
+  void     StrLinVertexFinderMinDist(Int_t OptUseWeights=0);
+  static void GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t m[][3],Double_t *d);
+  static void GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t *sigmasq,Double_t m[][3],Double_t *d);
+  static Double_t GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t *x0);
+  static Double_t GetDeterminant3X3(Double_t matr[][3]);
+ 
   void     VertexFitter();
 
   ClassDef(AliITSVertexerTracks,1) // 3D Vertexing with ITS tracks 

diff --git a/ITS/AliITStrackV2.cxx b/ITS/AliITStrackV2.cxx
index d08e56af33d04d38a3e915d174d302e31077254f..ae042babdd838b3ab2e4c087b1ab0415b047942f 100644 (file)
--- a/ITS/AliITStrackV2.cxx
+++ b/ITS/AliITStrackV2.cxx
@@ -13,17 +13,18 @@
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-//-------------------------------------------------------------------------
+///////////////////////////////////////////////////////////////////////////
 //                Implementation of the ITS track class
 //
 //          Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
 //     dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
-//-------------------------------------------------------------------------
+///////////////////////////////////////////////////////////////////////////
 #include <TMath.h>
 
 #include "AliCluster.h"
 #include "AliESDtrack.h"
 #include "AliITStrackV2.h"
+#include "AliITSStrLine.h"
 
 ClassImp(AliITStrackV2)
 
@@ -210,6 +211,24 @@ GetGlobalXYZat(Double_t xk, Double_t &x, Double_t &y, Double_t &z) const {
   return 1;
 }
 
+//_____________________________________________________________________________
+void AliITStrackV2::ApproximateHelixWithLine(Double_t xk, AliITSStrLine *line)
+{
+  //------------------------------------------------------------
+  // Approximate the track (helix) with a straight line tangent to the
+  // helix in the point defined by r (F. Prino, prino@to.infn.it)
+  //------------------------------------------------------------
+  Double_t mom[3];
+  Double_t azim = TMath::ASin(fP2)+fAlpha;
+  Double_t theta = TMath::Pi()/2. - TMath::ATan(fP3);
+  mom[0] = TMath::Sin(theta)*TMath::Cos(azim);
+  mom[1] = TMath::Sin(theta)*TMath::Sin(azim);
+  mom[2] = TMath::Cos(theta);
+  Double_t pos[3];
+  GetGlobalXYZat(xk,pos[0],pos[1],pos[2]);
+  line->SetP0(pos);
+  line->SetCd(mom);
+}
 //_____________________________________________________________________________
 Double_t AliITStrackV2::GetPredictedChi2(const AliCluster *c) const 
 {

diff --git a/ITS/AliITStrackV2.h b/ITS/AliITStrackV2.h
index 7f03f989ad6e03f8ac5f960be8d98a92750e27cc..8a766ec46b7af117be7dc3fba87d929f164975cb 100644 (file)
--- a/ITS/AliITStrackV2.h
+++ b/ITS/AliITStrackV2.h
@@ -30,6 +30,7 @@
 #include "AliITSrecoV2.h"
 
 class AliESDtrack;
+class AliITSStrLine;
 
 //_____________________________________________________________________________
 class AliITStrackV2 : public AliKalmanTrack {
@@ -76,6 +77,7 @@ public:
   void GetExternalCovariance(Double_t cov[15]) const ;
   Int_t GetClusterIndex(Int_t i) const {return fIndex[i];}
   Int_t GetGlobalXYZat(Double_t r,Double_t &x,Double_t &y,Double_t &z) const;
+  void ApproximateHelixWithLine(Double_t xk, AliITSStrLine *line);
   Double_t GetPredictedChi2(const AliCluster *cluster) const;
    Int_t Invariant() const;
  

diff --git a/ITS/ITSrecLinkDef.h b/ITS/ITSrecLinkDef.h
index fa63785fe0a7509314fde28028a4afe606501bd1..d866e3cba484102b73265684bdb4b34fbc3f5395 100644 (file)
--- a/ITS/ITSrecLinkDef.h
+++ b/ITS/ITSrecLinkDef.h
@@ -52,6 +52,7 @@
 #pragma link C++ class  AliITSVertexerPPZ+;
 #pragma link C++ class  AliITSVertexerZ+;
 #pragma link C++ class  AliITSStrLine+;
+#pragma link C++ class  AliITSSimpleVertex+;
 #pragma link C++ class  AliITSVertexerTracks+;
 
 // Classes for neural tracking
@@ -75,6 +76,7 @@
 #pragma link C++ class AliITSClusterFinderV2SSD+;
 
 //beam test classes
+#pragma link C++ class  AliITSspdTestBeam+;
 #pragma link C++ class AliITSBeamTestDig+;
 #pragma link C++ class AliITSBeamTestDigSPD+;
 #pragma link C++ class AliITSBeamTestDigSDD+;

diff --git a/ITS/ITSsimLinkDef.h b/ITS/ITSsimLinkDef.h
index 176e7e56be78c04611e1c03c830b1e10ab474c25..f8311f8ef05338980458147e814f398dd5a12aab 100644 (file)
--- a/ITS/ITSsimLinkDef.h
+++ b/ITS/ITSsimLinkDef.h
@@ -32,17 +32,6 @@
 //#pragma link C++ class  AliITSv11GeomCable+;
 //#pragma link C++ class  AliITSv11GeomcableFlat+;
 //#pragma link C++ class  AliITSv11GeomcableRound+;
-//
-#pragma link C++ class  AliITSspdTestBeam+;
-//#pragma link C++ class  AliITSTestBeamData-;
-//#pragma link C++ class  AliITSspdTestBeamHeader-;
-//#pragma link C++ class  AliITSspdTestBeamTail-;
-//#pragma link C++ class  AliITSspdTestBeamBurst-;
-//#pragma link C++ class  AliITSspdTestBeamData-;
-//
-//#pragma link C++ class  AliITSMixture+;
-//#pragma link C++ class  AliITSGeoCable+;
-//
 #pragma link C++ class  AliITShit+;
 #pragma link C++ class  AliITSGeant3Geometry+;
 // Standard ITS detector class initilizers

diff --git a/ITS/libITSrec.pkg b/ITS/libITSrec.pkg
index e39ed6d5a0b27e007e055de5f940e270f137f6b4..fa790749dbebe9a047f6e51388524de58d838938 100644 (file)
--- a/ITS/libITSrec.pkg
+++ b/ITS/libITSrec.pkg
@@ -21,6 +21,7 @@ SRCS =        AliITSDetTypeRec.cxx \
                AliITSVertexer.cxx \
                AliITSVertexerIons.cxx \
                AliITSVertexerPPZ.cxx \
+               AliITSSimpleVertex.cxx \
                AliITSVertexerTracks.cxx \
                AliITSVertexerZ.cxx \
                AliITSVertexerFast.cxx \
@@ -50,6 +51,7 @@ SRCS =        AliITSDetTypeRec.cxx \
                AliITSBeamTestDigSPD.cxx \
                AliITSBeamTestDigSSD.cxx \
                AliITSBeamTestDigitizer.cxx\
+                AliITSspdTestBeam.cxx \
 
 HDRS:=  $(SRCS:.cxx=.h)
 

diff --git a/ITS/libITSsim.pkg b/ITS/libITSsim.pkg
index 111e98fda8d2ced8159481c18098c364312fc08f..958a0518ec5908208114fbb0f6fa8e69f9408f5c 100644 (file)
--- a/ITS/libITSsim.pkg
+++ b/ITS/libITSsim.pkg
@@ -26,7 +26,6 @@ SRCS =        AliITS.cxx \
                AliITSsDigitize.cxx \
                AliITSDigitizer.cxx \
                AliITSFDigitizer.cxx \
-      AliITSspdTestBeam.cxx \
       AliITSTrigger.cxx