[u/mrichter/AliRoot.git] / ITS / AliITSVertexerIons.cxx

#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 "AliITSLoader.h"
#include "AliITSRecPoint.h"
#include "AliGenerator.h"
#include "AliMagF.h"

#include <TH1.h>
#include <TF1.h>
#include <TCanvas.h>
#include <AliESDVertex.h>
#include <TObjArray.h>
#include <TObject.h>
#include <AliITSVertexerPPZ.h>
//////////////////////////////////////////////////////////////////////
// AliITSVertexerIons is a class for full 3D primary vertex         //
// finding optimized for Ion-Ion interactions                       //
//                                                                  // 
//                                                                  //
//                                                                  //
//                                                                  //
// Written by Giuseppe Lo Re and Francesco Riggi                    //
// Giuseppe.Lore@ct.infn.it                                         //
// Franco.Riggi@ct.infn.it                                          //
//                                                                  //
// Release date: Mar 2004                                           //
//                                                                  //
//                                                                  //       
//////////////////////////////////////////////////////////////////////

ClassImp(AliITSVertexerIons)


//______________________________________________________________________
  AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer() {
  // Default Constructor

  fITS = 0;
  SetNpThreshold();
  SetMaxDeltaPhi();
  SetMaxDeltaZ();
}

//______________________________________________________________________
AliITSVertexerIons::AliITSVertexerIons(TString fn):AliITSVertexer(fn) {
  // Standard constructor
  
  fITS = 0;
  SetNpThreshold();
  SetMaxDeltaPhi();
  SetMaxDeltaZ();}


//______________________________________________________________________
AliITSVertexerIons::~AliITSVertexerIons() {
  // Default Destructor
  fITS = 0;
}

//______________________________________________________________________
AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){

  fCurrentVertex = 0;

  AliRunLoader *rl = AliRunLoader::GetRunLoader();
  AliITSLoader* itsloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
  if(!fITS) {
    fITS =(AliITS *)gAlice->GetDetector("ITS");
    if(!fITS) {
      Error("FindVertexForCurrentEvent","AliITS object was not found");
      return fCurrentVertex;
    }
  }
  fITS->SetTreeAddress();
  AliITSgeom *g2 = fITS->GetITSgeom();
  TClonesArray  *recpoints = fITS->RecPoints();
  AliITSRecPoint *pnt;
  TTree *tr =  itsloader->TreeR();

  Int_t npoints=0;
  Int_t nopoints1=40000;
  Int_t nopoints2=40000;
  Float_t l[3], 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=g2->GetStartSPD();i<=g2->GetLastSPD();i++) {
    fITS->ResetRecPoints();
    tr->GetEvent(i);
    npoints = recpoints->GetEntries();
    for (Int_t 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);
      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++;
      }
      if(i>=80 && TMath::Abs(p[2])<14.35) {
		y2[np2]=p[1];
		x2[np2]=p[0];
		z2[np2]=p[2];
		np2++;
      }
    }
  }

  if(np1<fNpThreshold) {
    Warning("FindVertexForCurrentEvent","AliITSVertexerIons finder is not reliable for low multiplicity events. Switching to AliITSVertexerPPZ with default parameters...\n");
    Warning("FindVertexForCurrentEvent","N rec points = %d - Threshold is %d",np1,fNpThreshold);
    AliITSVertexerPPZ *dovert = new AliITSVertexerPPZ("default");
    fCurrentVertex =dovert->FindVertexForCurrentEvent(rl->GetEventNumber());
    delete dovert;
    return fCurrentVertex;
  }

  if(!np1 || !np2) {
    Error("FindVertexForCurrentEvent","No points in the pixer layers");
    return fCurrentVertex;
  }  
  if(fDebug>0)   cout << "N. points layer 1 and 2 : " << np1 << " " << np2 << endl;

  Double_t asparx = (mxpiu-mxmeno)/(mxpiu+mxmeno);
  Double_t aspary = (mypiu-mymeno)/(mypiu+mymeno);

  Double_t x0 = 2.86*asparx;
  Double_t y0 = 2.86*aspary;

  if(fDebug>0)   cout << "Rough xy vertex = " << x0 << " " << y0 << endl;  

  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]);
  }
   
  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));
     }
  }

  // Fitting ...
  Double_t max;
  Int_t bin_max;
  Double_t max_center;
  
  max = hzv->GetMaximum();
  bin_max=hzv->GetMaximumBin();
  max_center=hzv->GetBinCenter(bin_max);
  Double_t dxy=1.5;
  
  TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",max_center-dxy,max_center+dxy);  
  fz->SetParameter(0,max);
  fz->SetParameter(1,max_center);
  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];
    if(fDebug>0)Info("FindVertexForCurrentEvent","Vertex found for event %d",evnumber);
    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 x_approx=FindMaxAround(x0,hxv,dist);
  if(fDebug>0) cout << "x_approx = " << x_approx << endl;
  fx->SetParameter(1,x_approx);
  Double_t dif_centroid=0.07;
  fx->SetParLimits(1,x_approx-dif_centroid,x_approx+dif_centroid);
  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 y_approx=FindMaxAround(y0,hyv,dist);
  if(fDebug>0) cout << "y_approx = " << y_approx << endl;
  fy->SetParameter(1,y_approx);
  fy->SetParLimits(1,y_approx-dif_centroid,y_approx+dif_centroid);
  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];
  if(fDebug>0)Info("FindVertexForCurrentEvent","Vertex found for event %d",evnumber);
  sprintf(name,"Vertex");
  fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
  return fCurrentVertex;
}

//______________________________________________________________________
void AliITSVertexerIons::PhiFunc(Double_t &x,Double_t &y,Double_t &phi) {
  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::FindVertices(){
  // computes the vertices of the events in the range FirstEvent - LastEvent
  AliRunLoader *rl = AliRunLoader::GetRunLoader();
  AliITSLoader* itsloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
  itsloader->LoadRecPoints("read");
  for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
    rl->GetEvent(i);
    FindVertexForCurrentEvent(i);
    if(fCurrentVertex){
      WriteCurrentVertex();
    }
    else {
      if(fDebug>0){
	cout<<"Vertex not found for event "<<i<<endl;
      }
    }
  }
}

//________________________________________________________
void AliITSVertexerIons::PrintStatus() const {
  // Print current status
  cout <<"=======================================================\n";
  cout <<" Debug flag: "<<fDebug<<endl;
  cout<<"First event to be processed "<<fFirstEvent;
  cout<<"\n Last event to be processed "<<fLastEvent<<endl;
  if(fCurrentVertex)fCurrentVertex->PrintStatus();
}

//________________________________________________________
Double_t AliITSVertexerIons::FindMaxAround(Double_t point, TH1F *h, Double_t distance) {
  Int_t max_content=0;
  Int_t max_bin=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(fabs(center-point)>distance) continue;
    if(content>max_content) {max_content=content;max_bin=i;}
  }
  Double_t max=h->GetBinCenter(max_bin);
  return max;
}
Commit	Line	Data
c5f0f3c1	1	#include <AliITSVertexerIons.h>
	2	#include "stdlib.h"
	3	#include <TMath.h>
	4	#include <TRandom.h>
	5	#include <TObjArray.h>
	6	#include <TROOT.h>
	7	#include <TFile.h>
	8	#include <TTree.h>
	9	#include <Riostream.h>
	10
	11	#include "AliRun.h"
	12	#include "AliITS.h"
	13	#include "AliITSgeom.h"
88cb7938	14	#include "AliITSLoader.h"
c5f0f3c1	15	#include "AliITSRecPoint.h"
	16	#include "AliGenerator.h"
	17	#include "AliMagF.h"
	18
	19	#include <TH1.h>
	20	#include <TF1.h>
	21	#include <TCanvas.h>
d681bb2d	22	#include <AliESDVertex.h>
c5f0f3c1	23	#include <TObjArray.h>
c5f0f3c1	24	#include <TObject.h>
83028144	25	#include <AliITSVertexerPPZ.h>
c5f0f3c1	26	//////////////////////////////////////////////////////////////////////
	27	// AliITSVertexerIons is a class for full 3D primary vertex //
	28	// finding optimized for Ion-Ion interactions //
	29	// //
	30	// //
	31	// //
	32	// //
	33	// Written by Giuseppe Lo Re and Francesco Riggi //
	34	// Giuseppe.Lore@ct.infn.it //
c5f0f3c1	35	// Franco.Riggi@ct.infn.it //
c5f0f3c1	36	// //
fe4280bf	37	// Release date: Mar 2004 //
c5f0f3c1	38	// //
	39	// //
	40	//////////////////////////////////////////////////////////////////////
	41
	42	ClassImp(AliITSVertexerIons)
	43
	44
	45
fe4280bf	46	//______________________________________________________________________
88cb7938	47	AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer() {
c5f0f3c1	48	// Default Constructor
	49
	50	fITS = 0;
83028144	51	SetNpThreshold();
fe4280bf	52	SetMaxDeltaPhi();
fe4280bf	53	SetMaxDeltaZ();
c5f0f3c1	54	}
	55
	56	//______________________________________________________________________
88cb7938	57	AliITSVertexerIons::AliITSVertexerIons(TString fn):AliITSVertexer(fn) {
c5f0f3c1	58	// Standard constructor
88cb7938	59
c5f0f3c1	60	fITS = 0;
83028144	61	SetNpThreshold();
fe4280bf	62	SetMaxDeltaPhi();
fe4280bf	63	SetMaxDeltaZ();}
c5f0f3c1	64
	65
	66	//______________________________________________________________________
	67	AliITSVertexerIons::~AliITSVertexerIons() {
	68	// Default Destructor
	69	fITS = 0;
	70	}
	71
	72	//______________________________________________________________________
d681bb2d	73	AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
fe4280bf	74
c5f0f3c1	75	fCurrentVertex = 0;
fe4280bf	76
88cb7938	77	AliRunLoader *rl = AliRunLoader::GetRunLoader();
83028144	78	AliITSLoader* itsloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
fe4280bf	79	if(!fITS) {
c5f0f3c1	80	fITS =(AliITS *)gAlice->GetDetector("ITS");
	81	if(!fITS) {
	82	Error("FindVertexForCurrentEvent","AliITS object was not found");
	83	return fCurrentVertex;
	84	}
	85	}
88cb7938	86	fITS->SetTreeAddress();
fe4280bf	87	AliITSgeom *g2 = fITS->GetITSgeom();
88cb7938	88	TClonesArray *recpoints = fITS->RecPoints();
88cb7938	89	AliITSRecPoint *pnt;
fe4280bf	90	TTree *tr = itsloader->TreeR();
88cb7938	91
fe4280bf	92	Int_t npoints=0;
	93	Int_t nopoints1=40000;
	94	Int_t nopoints2=40000;
88cb7938	95	Float_t l[3], p[3];
c5f0f3c1	96
83028144	97	Double_t z1, z2, y1, y2, x1, x2, phi1, phi2, r1, r2;
	98	z1=new Double_t[nopoints1];
	99	z2=new Double_t[nopoints2];
	100	y1=new Double_t[nopoints1];
	101	y2=new Double_t[nopoints2];
	102	x1=new Double_t[nopoints1];
	103	x2=new Double_t[nopoints2];
	104	phi1=new Double_t[nopoints1];
	105	phi2=new Double_t[nopoints2];
	106	r1=new Double_t[nopoints1];
	107	r2=new Double_t[nopoints2];
c5f0f3c1	108
fe4280bf	109	Double_t mxpiu = 0;
	110	Double_t mxmeno = 0;
	111	Double_t mypiu = 0;
	112	Double_t mymeno = 0;
	113	Double_t r=0;
	114
	115	Int_t np1=0, np2=0;
	116	for(Int_t i=g2->GetStartSPD();i<=g2->GetLastSPD();i++) {
	117	fITS->ResetRecPoints();
	118	tr->GetEvent(i);
	119	npoints = recpoints->GetEntries();
	120	for (Int_t ipoint=0;ipoint<npoints;ipoint++) {
	121	pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
	122	l[0]=pnt->GetX();
	123	l[1]=0;
	124	l[2]=pnt->GetZ();
	125	g2->LtoG(i, l, p);
	126	r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
	127
	128	if(i<80 && TMath::Abs(p[2])<14.35) {
	129	y1[np1]=p[1];
	130	x1[np1]=p[0];
	131	z1[np1]=p[2];
	132	if(p[0]>0) mxpiu++;
	133	if(p[0]<0) mxmeno++;
	134	if(p[1]>0) mypiu++;
	135	if(p[1]<0) mymeno++;
	136	np1++;
2257f27e	137	}
fe4280bf	138	if(i>=80 && TMath::Abs(p[2])<14.35) {
	139	y2[np2]=p[1];
	140	x2[np2]=p[0];
	141	z2[np2]=p[2];
	142	np2++;
88cb7938	143	}
88cb7938	144	}
fe4280bf	145	}
c5f0f3c1	146
83028144	147	if(np1<fNpThreshold) {
83028144	148	Warning("FindVertexForCurrentEvent","AliITSVertexerIons finder is not reliable for low multiplicity events. Switching to AliITSVertexerPPZ with default parameters...\n");
83028144	149	Warning("FindVertexForCurrentEvent","N rec points = %d - Threshold is %d",np1,fNpThreshold);
fe4280bf	150	AliITSVertexerPPZ *dovert = new AliITSVertexerPPZ("default");
83028144	151	fCurrentVertex =dovert->FindVertexForCurrentEvent(rl->GetEventNumber());
	152	delete dovert;
	153	return fCurrentVertex;
	154	}
	155
fe4280bf	156	if(!np1 \|\| !np2) {
fe4280bf	157	Error("FindVertexForCurrentEvent","No points in the pixer layers");
83028144	158	return fCurrentVertex;
fe4280bf	159	}
	160	if(fDebug>0) cout << "N. points layer 1 and 2 : " << np1 << " " << np2 << endl;
	161
	162	Double_t asparx = (mxpiu-mxmeno)/(mxpiu+mxmeno);
	163	Double_t aspary = (mypiu-mymeno)/(mypiu+mymeno);
	164
	165	Double_t x0 = 2.86*asparx;
	166	Double_t y0 = 2.86*aspary;
	167
	168	if(fDebug>0) cout << "Rough xy vertex = " << x0 << " " << y0 << endl;
	169
	170	for(Int_t i=0;i<np1;i++) {
	171	x1[i]-=x0;
	172	y1[i]-=y0;
	173	PhiFunc(x1[i],y1[i],phi1[i]);
	174	r1[i]=TMath::Sqrt(x1[i]x1[i]+y1[i]y1[i]);
88cb7938	175	}
fe4280bf	176	for(Int_t i=0;i<np2;i++) {
	177	x2[i]-=x0;
	178	y2[i]-=y0;
	179	PhiFunc(x2[i],y2[i],phi2[i]);
	180	r2[i]=TMath::Sqrt(x2[i]x2[i]+y2[i]y2[i]);
83028144	181	}
fe4280bf	182
	183	Int_t nbinxy=400;
	184	Int_t nbinz=400;
	185	TH1F *hxv=new TH1F("hxv","",nbinxy,-5,5);
	186	TH1F *hyv=new TH1F("hyv","",nbinxy,-5,5);
	187	TH1F *hzv=new TH1F("hzv","",nbinz,-15,15);
	188
	189	Double_t dphi;
	190	for(Int_t j=0; j<np1; j++) {
	191	for(Int_t k=0; k<np2; k++) {
	192	dphi=TMath::Abs(phi2[k]-phi1[j]);
88cb7938	193	if(dphi>180) dphi = 360-dphi;
fe4280bf	194	if(dphi>fMaxDeltaPhi) continue;
	195	hzv->Fill(z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1));
	196	}
88cb7938	197	}
fe4280bf	198
	199	// Fitting ...
	200	Double_t max;
	201	Int_t bin_max;
	202	Double_t max_center;
83028144	203
fe4280bf	204	max = hzv->GetMaximum();
	205	bin_max=hzv->GetMaximumBin();
	206	max_center=hzv->GetBinCenter(bin_max);
	207	Double_t dxy=1.5;
83028144	208
fe4280bf	209	TF1 fz = new TF1 ("fz","([0]exp(-0.5((x-[1])/[2])((x-[1])/[2])))+[3]",max_center-dxy,max_center+dxy);
	210	fz->SetParameter(0,max);
	211	fz->SetParameter(1,max_center);
	212	fz->SetParameter(2,0.1);
	213	fz->SetLineColor(kRed);
	214	hzv->Fit("fz","RQ0");
	215
	216	if(TMath::Abs(x0)>0.4 \|\| TMath::Abs(y0)>0.4) {
	217	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");
	218	Double_t position[3]={0,0,fz->GetParameter(1)};
	219	Double_t resolution[3]={0,0,0};
	220	Double_t snr[3]={0,0,0};
	221	Char_t name[30];
	222	if(fDebug>0)Info("FindVertexForCurrentEvent","Vertex found for event %d",evnumber);
	223	sprintf(name,"Vertex");
	224	fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
83028144	225	return fCurrentVertex;
88cb7938	226	}
83028144	227
fe4280bf	228	for(Int_t j=0; j<np1; j++) {
	229	for(Int_t k=0; k<np2; k++) {
	230	if(TMath::Abs((z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1))-fz->GetParameter(1))>fMaxDeltaZ) continue;
	231	if(y2[k]==y1[j]) continue;
	232	hxv->Fill(x0+(x2[k]-((x2[k]-x1[j])/(y2[k]-y1[j]))*y2[k]));
	233	if(x2[k]==x1[j]) continue;
	234	hyv->Fill(y0+(y2[k]-((y2[k]-y1[j])/(x2[k]-x1[j]))*x2[k]));
88cb7938	235	}
88cb7938	236	}
fe4280bf	237
	238	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);
	239	fx->SetParameter(0,100);
	240	Double_t dist=0.3;
	241	Double_t x_approx=FindMaxAround(x0,hxv,dist);
	242	if(fDebug>0) cout << "x_approx = " << x_approx << endl;
	243	fx->SetParameter(1,x_approx);
	244	Double_t dif_centroid=0.07;
	245	fx->SetParLimits(1,x_approx-dif_centroid,x_approx+dif_centroid);
	246	fx->SetParameter(2,0.1);
	247	fx->SetLineColor(kRed);
	248	hxv->Fit("fx","RQW0");
83028144	249
fe4280bf	250	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);
	251	fy->SetParameter(0,100);
	252	Double_t y_approx=FindMaxAround(y0,hyv,dist);
	253	if(fDebug>0) cout << "y_approx = " << y_approx << endl;
	254	fy->SetParameter(1,y_approx);
	255	fy->SetParLimits(1,y_approx-dif_centroid,y_approx+dif_centroid);
	256	fy->SetParameter(2,0.1);
	257	fy->SetLineColor(kRed);
	258	hyv->Fit("fy","RQW0");
83028144	259
83028144	260	delete [] z1;
	261	delete [] z2;
	262	delete [] y1;
	263	delete [] y2;
	264	delete [] x1;
	265	delete [] x2;
88cb7938	266	delete [] r1;
	267	delete [] r2;
	268	delete [] phi1;
	269	delete [] phi2;
fe4280bf	270	delete hxv;
	271	delete hyv;
	272	delete hzv;
	273
	274	Double_t position[3]={fx->GetParameter(1),fy->GetParameter(1),fz->GetParameter(1)};
	275	Double_t resolution[3]={fx->GetParameter(2),fy->GetParameter(2),fz->GetParameter(2)};
	276	Double_t snr[3]={0,0,0};
	277
88cb7938	278	Char_t name[30];
088e0b8d	279	if(fDebug>0)Info("FindVertexForCurrentEvent","Vertex found for event %d",evnumber);
88cb7938	280	sprintf(name,"Vertex");
d681bb2d	281	fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
88cb7938	282	return fCurrentVertex;
c5f0f3c1	283	}
	284
	285	//______________________________________________________________________
fe4280bf	286	void AliITSVertexerIons::PhiFunc(Double_t &x,Double_t &y,Double_t &phi) {
	287	if(y>0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578);
	288	if(y>0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180;
	289	if(y<0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180;
	290	if(y<0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+360;;
c5f0f3c1	291	}
	292
	293	//______________________________________________________________________
	294	void AliITSVertexerIons::FindVertices(){
	295	// computes the vertices of the events in the range FirstEvent - LastEvent
88cb7938	296	AliRunLoader *rl = AliRunLoader::GetRunLoader();
83028144	297	AliITSLoader* itsloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
83028144	298	itsloader->LoadRecPoints("read");
ab6a6f40	299	for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
88cb7938	300	rl->GetEvent(i);
c5f0f3c1	301	FindVertexForCurrentEvent(i);
	302	if(fCurrentVertex){
	303	WriteCurrentVertex();
	304	}
	305	else {
	306	if(fDebug>0){
	307	cout<<"Vertex not found for event "<<i<<endl;
	308	}
	309	}
	310	}
c5f0f3c1	311	}
	312
	313	//________________________________________________________
	314	void AliITSVertexerIons::PrintStatus() const {
	315	// Print current status
88cb7938	316	cout <<"=======================================================\n";
88cb7938	317	cout <<" Debug flag: "<<fDebug<<endl;
c5f0f3c1	318	cout<<"First event to be processed "<<fFirstEvent;
	319	cout<<"\n Last event to be processed "<<fLastEvent<<endl;
	320	if(fCurrentVertex)fCurrentVertex->PrintStatus();
	321	}
fe4280bf	322
	323	//________________________________________________________
	324	Double_t AliITSVertexerIons::FindMaxAround(Double_t point, TH1F *h, Double_t distance) {
	325	Int_t max_content=0;
	326	Int_t max_bin=0;
	327	for(Int_t i=0;i<h->GetNbinsX();i++) {
	328	Int_t content=(Int_t)h->GetBinContent(i);
	329	Double_t center=(Double_t)h->GetBinCenter(i);
	330	if(fabs(center-point)>distance) continue;
	331	if(content>max_content) {max_content=content;max_bin=i;}
	332	}
	333	Double_t max=h->GetBinCenter(max_bin);
	334	return max;
	335	}
	336
	337