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

#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 "AliLog.h"
#include <TTree.h>
#include <TH1.h>
#include <TF1.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(),
fNpThreshold(0),
fMaxDeltaPhi(0),
fMaxDeltaZ(0){
  // Default Constructor

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

//______________________________________________________________________
AliITSVertexerIons::AliITSVertexerIons(TString fn):AliITSVertexer(fn),
fNpThreshold(0),
fMaxDeltaPhi(0),
fMaxDeltaZ(0) {
  // Standard constructor
  
  //fITS = 0;
  SetNpThreshold();
  SetMaxDeltaPhi();
  SetMaxDeltaZ();
}
/*
//______________________________________________________________________
AliITSVertexerIons::AliITSVertexerIons(const AliITSVertexerIons &source):AliITSVertexer(source) {
  // Copy constructor
  // Copies are not allowed. The method is protected to avoid misuse.
  Error("AliITSVertexerIons","Copy constructor not allowed\n");
}
*/
//_________________________________________________________________________
//AliITSVertexerIons& AliITSVertexerIons::operator=(const AliITSVertexerIons &/*source*/) {
  // Assignment operator
  // Assignment is not allowed. The method is protected to avoid misuse.
  //Error("= operator","Assignment operator not allowed\n");
  //return *this;
//}


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

//______________________________________________________________________
AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){ 
// Defines the AliESDVertex for the current event

  fCurrentVertex = 0;

  AliRunLoader *rl = AliRunLoader::GetRunLoader();
  AliITSLoader* itsloader = (AliITSLoader*)rl->GetLoader("ITSLoader");
  /*
  TDirectory * olddir = gDirectory;
  rl->CdGAFile();
  AliITSgeom* g2  = (AliITSgeom*)gDirectory->Get("AliITSgeom");
  olddir->cd(); 
  */

  TTree *tr =  itsloader->TreeR();
  AliITSDetTypeRec detTypeRec;

  detTypeRec.SetTreeAddressR(tr);

  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();
    tr->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++;
      }
      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 AliITSVertexerZ with default parameters...\n");
    Warning("FindVertexForCurrentEvent","N rec points = %d - Threshold is %d",np1,fNpThreshold);
    AliITSVertexerZ *dovert = new AliITSVertexerZ("default");
    fCurrentVertex =dovert->FindVertexForCurrentEvent(rl->GetEventNumber());
    delete dovert;
    return fCurrentVertex;
  }

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

  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;

  AliDebug(1,Form("Rough xy vertex = %f %f",x0,y0));  

  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 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];
    AliDebug(1,Form("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 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];
  AliDebug(1,Form("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) {
  // 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::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 {
      AliDebug(1,Form("Vertex not found for event %d",i));
    }
  }
}

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