X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITSVertexerIons.cxx;h=65f3ddb0e967b4f84ca5c3fc7a465fa6dd7ee153;hb=b5a9f753d92077cec01150463deccec04904304a;hp=9738e7c0c5df0184e317b9ceb26a62d2d05bb90c;hpb=c5f0f3c15afed663b2477f5cebc9515fc56c727c;p=u%2Fmrichter%2FAliRoot.git

diff --git a/ITS/AliITSVertexerIons.cxx b/ITS/AliITSVertexerIons.cxx
index 9738e7c0c5d..65f3ddb0e96 100644
--- a/ITS/AliITSVertexerIons.cxx
+++ b/ITS/AliITSVertexerIons.cxx
@@ -1,26 +1,16 @@
-#include <AliITSVertexerIons.h>
-#include "stdlib.h"
-#include <TMath.h>
-#include <TRandom.h>
-#include <TObjArray.h>
-#include <TROOT.h>
-#include <TFile.h>
-#include <TTree.h>
-#include <Riostream.h>
-
-#include "AliRun.h"
-#include "AliITS.h"
-#include "AliITSgeom.h"
+#include "AliRunLoader.h"
+#include "AliITSDetTypeRec.h"
+#include "AliITSVertexerIons.h"
+#include "AliITSVertexerZ.h"
+#include "AliESDVertex.h"
+#include "AliITSgeomTGeo.h"
+#include "AliITSLoader.h"
 #include "AliITSRecPoint.h"
-#include "AliGenerator.h"
-#include "AliMagF.h"
-
+#include "AliLog.h"
+#include <TTree.h>
 #include <TH1.h>
 #include <TF1.h>
-#include <TCanvas.h>
-#include <AliITSVertex.h>
-#include <TObjArray.h>
-#include <TObject.h>
+
 //////////////////////////////////////////////////////////////////////
 // AliITSVertexerIons is a class for full 3D primary vertex         //
 // finding optimized for Ion-Ion interactions                       //
@@ -30,483 +20,273 @@
 //                                                                  //
 // Written by Giuseppe Lo Re and Francesco Riggi                    //
 // Giuseppe.Lore@ct.infn.it                                         //
-//                                                                  //
 // Franco.Riggi@ct.infn.it                                          //
 //                                                                  //
-// Release date: May 2001                                           //
+// Release date: Mar 2004                                           //
 //                                                                  //
 //                                                                  //       
 //////////////////////////////////////////////////////////////////////
 
-ClassImp(AliITSVertexerIons)
-
-
 
+ClassImp(AliITSVertexerIons)
+  
 //______________________________________________________________________
-AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer() {
+  AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer(),
+fNpThreshold(0),
+fMaxDeltaPhi(0),
+fMaxDeltaZ(0){
   // Default Constructor
 
-  fITS = 0;
+  //fITS = 0;
+  SetNpThreshold();
+  SetMaxDeltaPhi();
+  SetMaxDeltaZ();
 }
 
-//______________________________________________________________________
-AliITSVertexerIons::AliITSVertexerIons(TFile *infile, TFile *outfile):AliITSVertexer(infile,outfile) {
-  // Standard constructor
-  if(!fInFile)Fatal("AliITSVertexerIons","No inputfile has been defined!");
-  fITS = 0;
-}
-
-
 //______________________________________________________________________
 AliITSVertexerIons::~AliITSVertexerIons() {
   // Default Destructor
-  fITS = 0;
+  //fITS = 0;
 }
 
 //______________________________________________________________________
-AliITSVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
-  // Defines the AliITSVertex for the current event
-  fCurrentVertex = 0;
-  Double_t Position[3];
-  Double_t Resolution[3];
-  Double_t SNR[3];
- 
-  if(!fITS) { 
-    fITS =(AliITS *)gAlice->GetDetector("ITS");
-    if(!fITS) {
-      Error("FindVertexForCurrentEvent","AliITS object was not found");
-      return fCurrentVertex;
-    }
-  }
-   AliITSgeom *g2 = fITS->GetITSgeom(); 
-   
-   TClonesArray  *recpoints = fITS->RecPoints();
-   AliITSRecPoint *pnt;
+AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(TTree *itsClusterTree){ 
+// Defines the AliESDVertex for the current event
 
-   Int_t NoPoints1 = 0;                                         
-   Int_t NoPoints2 = 0;
-   Double_t Vzero[3];
-   Double_t AsPar[2];
-     
-//------------ Rough Vertex evaluation ---------------------------------   
-
-   Int_t i,npoints,ipoint,j,k,max,BinMax;
-   Double_t ZCentroid;
-   Float_t l[3], p[3];
-
-   Double_t mxpiu = 0;
-   Double_t mxmeno = 0;
-   Double_t mypiu = 0;
-   Double_t mymeno = 0;
-   
-   TH1F *hITSz1         = new TH1F("hITSz1","",100,-14.35,14.35);
-  
-   for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++) 
-   {
-           fITS->ResetRecPoints(); 
-       gAlice->TreeR()->GetEvent(i);    
-           npoints = recpoints->GetEntries();
-           
-           for (ipoint=0;ipoint<npoints;ipoint++) {
-               pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
-               l[0]=pnt->GetX();
-               l[1]=0;
-               l[2]=pnt->GetZ();
-                   g2->LtoG(i, l, p);
-                   if(i<80 && TMath::Abs(p[2])<14.35) {         
-                        if(p[0]>0) mxpiu++; 
-                        if(p[0]<0) mxmeno++; 
-                        if(p[1]>0) mypiu++; 
-                        if(p[1]<0) mymeno++; 
-                   hITSz1->Fill(p[2]);       
-               }
-                   if(i>=80 && TMath::Abs(p[2])<14.35) NoPoints2++;
-       }       
-   }  
-  
-   NoPoints1 = (Int_t)(hITSz1->GetEntries()); 
-           
-   AsPar[0] = (mxpiu-mxmeno)/(mxpiu+mxmeno);
-   AsPar[1] = (mypiu-mymeno)/(mypiu+mymeno); 
-   
-   Vzero[0] = 5.24441*AsPar[0]; 
-   Vzero[1] = 5.24441*AsPar[1]; 
-
-   ZCentroid=  TMath::Abs(hITSz1->GetMean()); 
-   Vzero[2] = -5.31040e-02+1.42198*ZCentroid+7.44718e-01*TMath::Power(ZCentroid,2)
-                  -5.73426e-01*TMath::Power(ZCentroid,3)+2.01500e-01*TMath::Power(ZCentroid,4)          
-              -3.34118e-02*TMath::Power(ZCentroid,5)+2.20816e-03*TMath::Power(ZCentroid,6);
-   
-   if(hITSz1->GetMean()<0) Vzero[2] = -Vzero[2];
-   
-   /*cout << "\nXvzero: " << Vzero[0] << " cm" << "";
-   cout << "\nYvzero: " << Vzero[1] << " cm" << "";
-   cout << "\nZvzero: " << Vzero[2] << " cm" << "\n";*/
-
-   delete hITSz1;
-
-   Double_t dphi,r,DeltaZ,a,b;
-   Int_t np1=0;
-   Int_t np2=0;
-   Int_t niter=0;
-   
-   Double_t DeltaPhiZ = 0.08;
-   
-   Float_t B[3];
-   Float_t origin[3];
-   for(Int_t okm=0;okm<3;okm++) origin[okm]=0;
-   gAlice->Field()->Field(origin,B);
-
-   DeltaPhiZ = DeltaPhiZ*B[2]/2;
-   Double_t DeltaPhiXY = 1.0;   
-
-//   cout << "\nDeltaPhiZ: " << DeltaPhiZ << " deg" << "\n";   
-//   cout << "DeltaPhiXY: " << DeltaPhiXY << " deg" << "\n";   
-   
-   Double_t *Z1, *Z2, *Y1, *Y2, *X1, *X2, *phi1, *phi2, *r1, *r2;
-   Z1=new Double_t[NoPoints1];
-   Z2=new Double_t[NoPoints2];
-   Y1=new Double_t[NoPoints1];
-   Y2=new Double_t[NoPoints2];
-   X1=new Double_t[NoPoints1];
-   X2=new Double_t[NoPoints2];
-   phi1=new Double_t[NoPoints1];
-   phi2=new Double_t[NoPoints2];
-   r1=new Double_t[NoPoints1];
-   r2=new Double_t[NoPoints2];
-        
-   DeltaZ = 4.91617e-01+2.67567e-02*Vzero[2]+1.49626e-02*TMath::Power(Vzero[2],2); 
-   Float_t MulFactorZ = 28000./(Float_t)NoPoints1;
-   Int_t nbin=(Int_t)((DeltaZ/0.005)/MulFactorZ);
-   Int_t nbinxy=250;
-   Int_t *VectorBinZ,*VectorBinXY;
-   VectorBinZ=new Int_t[nbin];
-   VectorBinXY=new Int_t[nbinxy];
-   Float_t f1= 0;
-   Float_t f2= 0;
-   Double_t sigma,MediaFondo;
-     
-   TH1D *hITSZv         = new TH1D("hITSZv","",nbin,Vzero[2]-DeltaZ,Vzero[2]+DeltaZ);
-   TH1D *hITSXv         = new TH1D("hITSXv","",nbinxy,-3,3);
-   TH1D *hITSYv         = new TH1D("hITSYv","",nbinxy,-3,3);
-
-//   cout << "DeltaZeta: " << DeltaZ << " cm" << "\n";   
-   
-   
-   start:   
-
-   hITSZv->Add(hITSZv,-1.);
-   hITSXv->Add(hITSXv,-1.);
-   hITSYv->Add(hITSYv,-1.);
-
-   np1=np2=0;
-
-   
-
-   for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++) 
-   {
-        fITS->ResetRecPoints(); 
-        gAlice->TreeR()->GetEvent(i);
-           npoints = recpoints->GetEntries();
-           for (ipoint=0;ipoint<npoints;ipoint++) {
-               
-                    pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
-                l[0]=pnt->GetX();
-                l[1]=0;
-                l[2]=pnt->GetZ();
-                g2->LtoG(i, l, p);
-              
-                        if(i<80 && TMath::Abs(p[2])<14.35) {
-                    p[0]=p[0]-Vzero[0];
-                    p[1]=p[1]-Vzero[1];
-                r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
-                        Y1[np1]=p[1];
-                        X1[np1]=p[0];
-                        Z1[np1]=p[2]; 
-                        r1[np1]=r;
-                        phi1[np1]=PhiFunc(p);
-                        np1++;
-                        }
-
-                        if(i>=80 &&  TMath::Abs(p[2])<14.35) { 
-                    p[0]=p[0]-Vzero[0];
-                    p[1]=p[1]-Vzero[1];
-                r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
-                        Y2[np2]=p[1];
-                        X2[np2]=p[0];
-                        Z2[np2]=p[2];
-                        r2[np2]=r;
-                        phi2[np2]=PhiFunc(p);
-                        np2++;
-                        }
-                        
-          }
-   }
-
-//------------------ Correlation between rec points ----------------------
-    
-   //cout << "\nNo. of Points on the two pixel layers: "<<np1<<" "<<np2<<endl; 
+  fCurrentVertex = 0;
 
-   for(j=0; j<(np2)-1; j++) {
-         for(k=0; k<(np1)-1; k++) { 
-                dphi=TMath::Abs(phi1[k]-phi2[j]);
-                if(dphi>180) dphi = 360-dphi;
-                if(dphi<DeltaPhiZ && TMath::Abs((Z1[k]-(Z2[j]-Z1[k])/((r2[j]/r1[k])-1))-Vzero[2])
-                         <DeltaZ) hITSZv->Fill(Z1[k]-(Z2[j]-Z1[k])/((r2[j]/r1[k])-1));        
-         }
-   }
-      
-   //cout << "\nNumber of used pairs: \n";
-   //cout << hITSZv->GetEntries() << '\n' << '\n';      
-   a = Vzero[2]-DeltaZ; 
-   b = Vzero[2]+DeltaZ; 
-   max=(Int_t) hITSZv->GetMaximum();
-   BinMax=hITSZv->GetMaximumBin();
-   sigma=0;
-   for(i=0;i<nbin;i++) VectorBinZ[i]=(Int_t)hITSZv->GetBinContent(i);
-   for(i=0;i<10;i++) f1=f1+VectorBinZ[i]/10;
-   for(i=nbin-10;i<nbin;i++) f2=f2+VectorBinZ[i]/10;
-   MediaFondo=(f1+f2)/2;
-   for(i=0;i<nbin;i++) {
-	  if(VectorBinZ[i]-MediaFondo>(max-MediaFondo)*0.4 && 
-	  VectorBinZ[i]-MediaFondo<(max-MediaFondo)*0.7) {
-		 sigma=hITSZv->GetBinCenter(BinMax)-hITSZv->GetBinCenter(i);
-		 sigma=TMath::Abs(sigma);
-		 if(sigma==0) sigma=0.05;
+  AliITSDetTypeRec detTypeRec;
+
+  detTypeRec.SetTreeAddressR(itsClusterTree);
+
+  TClonesArray  *recpoints = detTypeRec.RecPoints();
+  AliITSRecPoint *pnt;
+
+
+  Int_t npoints=0;
+  Int_t nopoints1=40000;
+  Int_t nopoints2=40000;
+  Float_t p[3];
+
+  Double_t *z1, *z2, *y1, *y2, *x1, *x2, *phi1, *phi2, *r1, *r2;
+  z1=new Double_t[nopoints1];
+  z2=new Double_t[nopoints2];
+  y1=new Double_t[nopoints1];
+  y2=new Double_t[nopoints2];
+  x1=new Double_t[nopoints1];
+  x2=new Double_t[nopoints2];
+  phi1=new Double_t[nopoints1];
+  phi2=new Double_t[nopoints2];
+  r1=new Double_t[nopoints1];
+  r2=new Double_t[nopoints2];
+
+  Double_t mxpiu = 0;
+  Double_t mxmeno = 0;
+  Double_t mypiu = 0;
+  Double_t mymeno = 0;
+  Double_t r=0;
+
+  Int_t np1=0, np2=0;
+  for(Int_t i=AliITSgeomTGeo::GetModuleIndex(1,1,1);i<=AliITSgeomTGeo::GetModuleIndex(2,1,1)-1;i++) {
+    detTypeRec.ResetRecPoints();
+    itsClusterTree->GetEvent(i);
+    npoints = recpoints->GetEntries();
+    for (Int_t ipoint=0;ipoint<npoints;ipoint++) {
+      pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
+      /*
+      l[0]=pnt->GetDetLocalX();
+      l[1]=0;
+      l[2]=pnt->GetDetLocalZ();
+      g2->LtoG(i, l, p);
+      */
+      pnt->GetGlobalXYZ(p);
+      r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
+
+      if(i<80 && TMath::Abs(p[2])<14.35)  {
+		y1[np1]=p[1];
+		x1[np1]=p[0];
+		z1[np1]=p[2];
+		if(p[0]>0) mxpiu++;
+		if(p[0]<0) mxmeno++;
+		if(p[1]>0) mypiu++;
+		if(p[1]<0) mymeno++;
+		np1++;
       }
-   }
-   
-   /*cout << "f1 " <<f1 <<endl;
-   cout << "f2 " <<f2 <<endl;
-   cout << "GetMaximumBin " <<hITSZv->GetMaximumBin() <<endl;
-   cout << "nbin " <<nbin <<endl;
-   cout << "max " << hITSZv->GetBinContent(BinMax)<<endl;
-   cout << "sigma " <<sigma <<endl;
-   cout << "Fondo " <<   MediaFondo<< endl;*/
-   
-   TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",a,b);  
-   
-   fz->SetParameter(0,max);
-   if(niter==0) {Double_t temp = hITSZv->GetBinCenter(BinMax); Vzero[2]=temp;}
-   fz->SetParameter(1,Vzero[2]);
-   fz->SetParameter(2,sigma); 
-   fz->SetParameter(3,MediaFondo);  
-   fz->SetParLimits(2,0,999);
-   fz->SetParLimits(3,0,999);
-   
-   hITSZv->Fit("fz","RMEQ0");
-     
-   SNR[2] = fz->GetParameter(0)/fz->GetParameter(3);
-   if(SNR[2]<0.) { 
-     cout << "\nNegative Signal to noise ratio for z!!!" << endl;
-     cout << "The algorithm cannot find the z vertex position." << endl;
-     exit(123456789);
-   }
-   else
-   {
-     Position[2] = fz->GetParameter(1);
-     if(Position[2]<0) 
-       {
-         Position[2]=Position[2]-TMath::Abs(Position[2])*1.11/10000;
-       }
-     else
-       {
-         Position[2]=Position[2]+TMath::Abs(Position[2])*1.11/10000;
-       }
-   }
-   Resolution[2] = fz->GetParError(1);
-   
-   
-   
-   for(j=0; j<(np2)-1; j++) {
-         for(k=0; k<(np1)-1; k++) { 
-                dphi=TMath::Abs(phi1[k]-phi2[j]);
-                if(dphi>180) dphi = 360-dphi;
-                if(dphi>DeltaPhiXY) continue;
-            if(TMath::Abs((Z1[k]-(Z2[j]-Z1[k])/((r2[j]/r1[k])-1))-Position[2])
-                <4*Resolution[2]) {
-                   hITSXv->Fill(Vzero[0]+(X2[j]-((X2[j]-X1[k])/(Y2[j]-Y1[k]))*Y2[j]));
-                   hITSYv->Fill(Vzero[1]+(Y2[j]-((Y2[j]-Y1[k])/(X2[j]-X1[k]))*X2[j]));
-            }
-         }
-   }
-   
-   TF1 *fx = new TF1
-   ("fx","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",Vzero[0]-0.5,Vzero[0]+0.5);  
-
-   max=(Int_t) hITSXv->GetMaximum();
-   BinMax=hITSXv->GetMaximumBin();
-   sigma=0;
-   f1=f2=0;
-   for(i=0;i<nbinxy;i++) VectorBinXY[i]=(Int_t)hITSXv->GetBinContent(i);
-   for(i=0;i<10;i++) f1=f1+VectorBinXY[i]/10;
-   for(i=nbinxy-10;i<nbinxy;i++) f2=f2+VectorBinXY[i]/10;
-   MediaFondo=(f1+f2)/2;
-   for(i=0;i<nbinxy;i++) {
-	  if(VectorBinXY[i]-MediaFondo>(max-MediaFondo)*0.4 && 
-	  VectorBinXY[i]-MediaFondo<(max-MediaFondo)*0.7) {
-		 sigma=hITSXv->GetBinCenter(BinMax)-hITSXv->GetBinCenter(i);
-		 sigma=TMath::Abs(sigma);
-		 if(sigma==0) sigma=0.05;
+      if(i>=80 && TMath::Abs(p[2])<14.35) {
+		y2[np2]=p[1];
+		x2[np2]=p[0];
+		z2[np2]=p[2];
+		np2++;
       }
-   }
-   
-   /*cout << "f1 " <<f1 <<endl;
-   cout << "f2 " <<f2 <<endl;
-   cout << "GetMaximumBin " <<hITSXv->GetMaximumBin() <<endl;
-   cout << "max " << hITSXv->GetBinContent(BinMax)<<endl;
-   cout << "sigma " <<sigma <<endl;
-   cout << "Fondo " <<   MediaFondo<< endl;
-   cout << "nbinxy " <<nbinxy <<endl;*/
-   
-   fx->SetParameter(0,max);
-   fx->SetParameter(1,Vzero[0]);
-   fx->SetParameter(2,sigma); 
-   fx->SetParameter(3,MediaFondo);
-
-   hITSXv->Fit("fx","RMEQ0"); 
-
-   SNR[0] = fx->GetParameter(0)/fx->GetParameter(3);
-   if(SNR[0]<0.) { 
-     cout << "\nNegative Signal to noise ratio for x!!!" << endl;
-     cout << "The algorithm cannot find the x vertex position." << endl;
-     exit(123456789);  
-   }
-   else
-   {
-   Position[0]=fx->GetParameter(1);
-   Resolution[0]=fx->GetParError(1);
-   }
+    }
+  }
 
-   TF1 *fy = new TF1("fy","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",Vzero[1]-0.5,Vzero[1]+0.5);  
+  if(np1<fNpThreshold) {
+    Warning("FindVertexForCurrentEvent","AliITSVertexerIons finder is not reliable for low multiplicity events. Switching to AliITSVertexerZ with default parameters...\n");
+    Warning("FindVertexForCurrentEvent","N rec points = %d - Threshold is %d",np1,fNpThreshold);
+    AliITSVertexerZ *dovert = new AliITSVertexerZ();
+    fCurrentVertex =dovert->FindVertexForCurrentEvent(itsClusterTree);
+    delete dovert;
+    return fCurrentVertex;
+  }
 
-   max=(Int_t) hITSYv->GetMaximum();
-   BinMax=hITSYv->GetMaximumBin();
-   sigma=0;
-   f1=f2=0;
-   for(i=0;i<nbinxy;i++) VectorBinXY[i]=(Int_t)hITSYv->GetBinContent(i);
-   for(i=0;i<10;i++) f1=f1+VectorBinXY[i]/10;
-   for(i=nbinxy-10;i<nbinxy;i++) f2=f2+VectorBinXY[i]/10;
-   MediaFondo=(f1+f2)/2;
-   for(i=0;i<nbinxy;i++) {
-	  if(VectorBinXY[i]-MediaFondo>(max-MediaFondo)*0.4 && 
-	  VectorBinXY[i]-MediaFondo<(max-MediaFondo)*0.7) {
-		 sigma=hITSYv->GetBinCenter(BinMax)-hITSYv->GetBinCenter(i);
-		 sigma=TMath::Abs(sigma);
-		 if(sigma==0) sigma=0.05;
-      }
-   }
-   
-   /*cout << "f1 " <<f1 <<endl;
-   cout << "f2 " <<f2 <<endl;
-   cout << "GetMaximumBin " <<hITSYv->GetMaximumBin() <<endl;
-   cout << "max " << hITSYv->GetBinContent(BinMax)<<endl;
-   cout << "sigma " <<sigma <<endl;
-   cout << "Fondo " <<   MediaFondo<< endl;
-   cout << "nbinxy " <<nbinxy <<endl;*/
-   
-   fy->SetParameter(0,max);
-   fy->SetParameter(1,Vzero[1]);
-   fy->SetParameter(2,sigma); 
-   fy->SetParameter(3,MediaFondo);
-   
-   hITSYv->Fit("fy","RMEQ0"); 
+  if(!np1 || !np2) {
+    Error("FindVertexForCurrentEvent","No points in the pixer layers");
+    return fCurrentVertex;
+  }  
+  AliDebug(1,Form("N. points layer 1 and 2 : %d %d",np1,np2));
 
-   SNR[1] = fy->GetParameter(0)/fy->GetParameter(3);
-   if(SNR[1]<0.) { 
-     cout << "\nNegative Signal to noise ratio for y!!!" << endl;
-     cout << "The algorithm cannot find the y vertex position." << endl;
-     exit(123456789); 
-   }
-   else
-   {
-   Position[1]=fy->GetParameter(1);
-   Resolution[1]=fy->GetParError(1);
-   }
+  Double_t asparx = (mxpiu-mxmeno)/(mxpiu+mxmeno);
+  Double_t aspary = (mypiu-mymeno)/(mypiu+mymeno);
 
-   niter++;
-   
-   Vzero[0] = Position[0];
-   Vzero[1] = Position[1];
-   Vzero[2] = Position[2];
+  Double_t x0 = 2.86*asparx;
+  Double_t y0 = 2.86*aspary;
 
-   if(niter<3) goto start;  // number of iterations
-   
+  AliDebug(1,Form("Rough xy vertex = %f %f",x0,y0));  
 
-   cout<<"Number of iterations: "<<niter<<endl;
- 
-   delete [] Z1;
-   delete [] Z2;
-   delete [] Y1;
-   delete [] Y2;
-   delete [] X1;
-   delete [] X2;
-   delete [] r1;
-   delete [] r2;
-   delete [] phi1;
-   delete [] phi2;
-   delete [] VectorBinZ;
-   delete [] VectorBinXY;
+  for(Int_t i=0;i<np1;i++) {
+    x1[i]-=x0;
+    y1[i]-=y0;
+    PhiFunc(x1[i],y1[i],phi1[i]);
+    r1[i]=TMath::Sqrt(x1[i]*x1[i]+y1[i]*y1[i]);
+  }
+  for(Int_t i=0;i<np2;i++) {
+    x2[i]-=x0;
+    y2[i]-=y0;
+    PhiFunc(x2[i],y2[i],phi2[i]);
+    r2[i]=TMath::Sqrt(x2[i]*x2[i]+y2[i]*y2[i]);
+  }
    
-   delete hITSZv;
-   delete hITSXv;
-   delete hITSYv;
-   Char_t name[30];
-   sprintf(name,"Vertex_%d",evnumber);
-   fCurrentVertex = new AliITSVertex(Position,Resolution,SNR,name);
-   return fCurrentVertex;
-}
+  Int_t nbinxy=400;
+  Int_t nbinz=400;
+  TH1F *hxv=new TH1F("hxv","",nbinxy,-5,5);
+  TH1F *hyv=new TH1F("hyv","",nbinxy,-5,5);
+  TH1F *hzv=new TH1F("hzv","",nbinz,-15,15);
+
+  Double_t dphi;
+  for(Int_t j=0; j<np1; j++) {
+    for(Int_t k=0; k<np2; k++) {
+      dphi=TMath::Abs(phi2[k]-phi1[j]);
+      if(dphi>180) dphi = 360-dphi;
+      if(dphi>fMaxDeltaPhi) continue; 
+      hzv->Fill(z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1));
+     }
+  }
 
-//______________________________________________________________________
-Double_t AliITSVertexerIons::PhiFunc(Float_t p[]) {
-    Double_t phi=0;
+  // Fitting ...
+  Double_t max;
+  Int_t binMax;
+  Double_t maxcenter;
+  
+  max = hzv->GetMaximum();
+  binMax=hzv->GetMaximumBin();
+  maxcenter=hzv->GetBinCenter(binMax);
+  Double_t dxy=1.5;
+  
+  TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",maxcenter-dxy,maxcenter+dxy);  
+  fz->SetParameter(0,max);
+  fz->SetParameter(1,maxcenter);
+  fz->SetParameter(2,0.1);
+  fz->SetLineColor(kRed);
+  hzv->Fit("fz","RQ0");
+
+  if(TMath::Abs(x0)>0.4 || TMath::Abs(y0)>0.4) {
+    Warning("FindVertexForCurrentEvent","AliITSVertexerIonsGeneral finder is not reliable for events with x and y vertex coordinates too far from the origin of the reference system. Their values will be set to 0.\n");
+    Double_t position[3]={0,0,fz->GetParameter(1)};
+    Double_t resolution[3]={0,0,0};
+    Double_t snr[3]={0,0,0};
+    Char_t name[30];
+    sprintf(name,"Vertex");
+    fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
+    return fCurrentVertex;
+  }
+  
+  for(Int_t j=0; j<np1; j++) {
+    for(Int_t k=0; k<np2; k++) {
+      if(TMath::Abs((z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1))-fz->GetParameter(1))>fMaxDeltaZ) continue;
+      if(y2[k]==y1[j]) continue;
+      hxv->Fill(x0+(x2[k]-((x2[k]-x1[j])/(y2[k]-y1[j]))*y2[k]));
+      if(x2[k]==x1[j]) continue;
+      hyv->Fill(y0+(y2[k]-((y2[k]-y1[j])/(x2[k]-x1[j]))*x2[k]));
+    }
+  }
+    
+  TF1 *fx = new TF1 ("fx","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]+[4]*x+[5]*x*x",x0-dxy,x0+dxy);  
+  fx->SetParameter(0,100);
+  Double_t dist=0.3;
+  Double_t xapprox=FindMaxAround(x0,hxv,dist);
+  AliDebug(1,Form("xapprox = %f",xapprox));
+  fx->SetParameter(1,xapprox);
+  Double_t difcentroid=0.07;
+  fx->SetParLimits(1,xapprox-difcentroid,xapprox+difcentroid);
+  fx->SetParameter(2,0.1);
+  fx->SetLineColor(kRed);
+  hxv->Fit("fx","RQW0"); 
+  
+  TF1 *fy = new TF1 ("fy","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]+[4]*x+[5]*x*x",y0-dxy,y0+dxy);  
+  fy->SetParameter(0,100);
+  Double_t yapprox=FindMaxAround(y0,hyv,dist);
+  AliDebug(1,Form("yapprox = %f",yapprox));
+  fy->SetParameter(1,yapprox);
+  fy->SetParLimits(1,yapprox-difcentroid,yapprox+difcentroid);
+  fy->SetParameter(2,0.1);
+  fy->SetLineColor(kRed);
+  hyv->Fit("fy","RQW0");
+  
+  delete [] z1;
+  delete [] z2;
+  delete [] y1;
+  delete [] y2;
+  delete [] x1;
+  delete [] x2;
+  delete [] r1;
+  delete [] r2;
+  delete [] phi1;
+  delete [] phi2;
+  delete hxv;
+  delete hyv;
+  delete hzv;
+  
+  Double_t position[3]={fx->GetParameter(1),fy->GetParameter(1),fz->GetParameter(1)};
+  Double_t resolution[3]={fx->GetParameter(2),fy->GetParameter(2),fz->GetParameter(2)};
+  Double_t snr[3]={0,0,0};
+  
+  Char_t name[30];
+  sprintf(name,"Vertex");
+  fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
 
-    if(p[1]>0 && p[0]>0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578);
-    if(p[1]>0 && p[0]<0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+180;
-    if(p[1]<0 && p[0]<0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+180;
-    if(p[1]<0 && p[0]>0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+360;
-    return phi;
+  return fCurrentVertex;
 }
 
 //______________________________________________________________________
-void AliITSVertexerIons::FindVertices(){
-  // computes the vertices of the events in the range FirstEvent - LastEvent
-  if(!fOutFile){
-    Error("FindVertices","The output file is not available - nothing done");
-    return;
-  }
-  if(!gAlice)gAlice = (AliRun*)fInFile->Get("gAlice");
-  if(!gAlice){
-    Error("FindVertices","The AliRun object is not available - nothing done");
-    return;
-  }
-  TDirectory *curdir = gDirectory;
-  fInFile->cd();
-  for(Int_t i=fFirstEvent;i<fLastEvent;i++){
-    gAlice->GetEvent(i);
-    FindVertexForCurrentEvent(i);
-    if(fCurrentVertex){
-      WriteCurrentVertex();
-    }
-    else {
-      if(fDebug>0){
-	cout<<"Vertex not found for event "<<i<<endl;
-      }
-    }
-  }
-  curdir->cd();
+void AliITSVertexerIons::PhiFunc(Double_t &x,Double_t &y,Double_t &phi) {
+  // Method for the computation of the phi angle given x and y. The result is in degrees. 
+  if(y>0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578);
+  if(y>0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180;
+  if(y<0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180;
+  if(y<0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+360;;
 }
 
 //________________________________________________________
 void AliITSVertexerIons::PrintStatus() const {
   // Print current status
- cout <<"=======================================================\n";
- cout <<" Debug flag: "<<fDebug<<endl;
-  if(fInFile)cout <<" The input file is open\n";
-  if(!fInFile)cout <<"The input file is not open\n";
-  if(fOutFile)cout <<"The output file is open\n";
-  if(fOutFile)cout <<"The output file not is open\n";
-  cout<<"First event to be processed "<<fFirstEvent;
-  cout<<"\n Last event to be processed "<<fLastEvent<<endl;
+  cout <<"=======================================================\n";
   if(fCurrentVertex)fCurrentVertex->PrintStatus();
 }
+
+//________________________________________________________
+Double_t AliITSVertexerIons::FindMaxAround(Double_t point, TH1F *h, Double_t distance) {
+  // It finds the maximum of h within point-distance and distance+point. 
+  Int_t maxcontent=0;
+  Int_t maxbin=0;
+  for(Int_t i=0;i<h->GetNbinsX();i++) {
+    Int_t content=(Int_t)h->GetBinContent(i);
+    Double_t center=(Double_t)h->GetBinCenter(i);
+    if(TMath::Abs(center-point)>distance) continue;
+    if(content>maxcontent) {maxcontent=content;maxbin=i;}
+  }
+  Double_t max=h->GetBinCenter(maxbin);
+  return max;
+}
+
+