[u/mrichter/AliRoot.git] / HLT / ITS / AliL3ITSVertexerZ.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, 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.                  *
 **************************************************************************/
#include "AliL3ITSVertexerZ.h"
#include<TString.h>
#include<TH1.h>
#include<TMath.h>
#include <AliRun.h>
#include <AliITS.h>
#include "AliITSLoader.h"
#include <AliITSgeom.h>
#include <AliITSRecPoint.h>

//-------------------------------------------------------------------------
//                Implementation of the HLT ITS vertexer class
//
//          Origin: Cvetan Cheshkov, CERN, Cvetan.Cheshkov@cern.ch
//-------------------------------------------------------------------------

ClassImp(AliL3ITSVertexerZ)

AliL3ITSVertexerZ::AliL3ITSVertexerZ(TString filename,Float_t x0, Float_t y0):AliITSVertexerZ(filename,x0,y0)
{
  // Standard Constructor
}

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

  fCurrentVertex = 0;
  AliRunLoader *rl =AliRunLoader::GetRunLoader();
  AliITSLoader* itsLoader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
  itsLoader->LoadRecPoints();
  rl->GetEvent(evnumber);

  if(!fITS)  {
    fITS = (AliITS*)gAlice->GetModule("ITS");
    if(!fITS) {
      Error("FindVertexForCurrentEvent","AliITS object was not found");
      return fCurrentVertex;
    }
  }

  fITS->SetTreeAddress();

  AliITSgeom *geom = fITS->GetITSgeom();

  TTree *tR = itsLoader->TreeR();
  TClonesArray *itsRec  = 0;
  Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
  Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
  Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
  Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;

  itsRec = fITS->RecPoints();

  //cout<<"Address of itsRec = "<<itsRec<<endl;
  TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
  TBranch *branch;
  if(fUseV2Clusters){
    branch = tR->GetBranch("Clusters");
    branch->SetAddress(&clusters);
  }
  else {
    branch = tR->GetBranch("ITSRecPoints");
  }


  Int_t nrpL1 = 0;
  Int_t nrpL2 = 0;
  for(Int_t module= fFirstL1; module<=fLastL1;module++){
    if(module%4==0 || module%4==3)continue;
    //   cout<<"Procesing module "<<module<<" ";
    branch->GetEvent(module);
    //    cout<<"Number of clusters "<<clusters->GetEntries()<<endl;
    if(fUseV2Clusters){
      Clusters2RecPoints(clusters,module,itsRec);
    }
    nrpL1+= itsRec->GetEntries();
    fITS->ResetRecPoints();
  }
  for(Int_t module= fFirstL2; module<=fLastL2;module++){
    branch->GetEvent(module);
    if(fUseV2Clusters){
      Clusters2RecPoints(clusters,module,itsRec);
    }
    nrpL2+= itsRec->GetEntries();
    fITS->ResetRecPoints();
  }
  if(nrpL1 == 0 || nrpL2 == 0){
    ResetHistograms();
    return fCurrentVertex;
  }

  Int_t nPhiBins = (Int_t)(TMath::Pi()/fDiffPhiMax);
  Float_t phiBinSize = 2*TMath::Pi()/(Float_t)nPhiBins;

  Int_t maxind1 = 2*nrpL1/nPhiBins;
  Float_t **zc1 = new Float_t *[nPhiBins];
  Float_t **phi1 = new Float_t *[nPhiBins];
  Float_t **r1 = new Float_t *[nPhiBins];
  Int_t *ind1 = new Int_t [nPhiBins];
  Int_t maxind2 = 2*nrpL2/nPhiBins;
  Float_t **zc2 = new Float_t *[nPhiBins];
  Float_t **phi2 = new Float_t *[nPhiBins];
  Float_t **r2 = new Float_t *[nPhiBins];
  Int_t *ind2 = new Int_t [nPhiBins];
  for(Int_t i=0;i<nPhiBins;i++) {
    zc1[i] = new Float_t [maxind1];
    phi1[i] = new Float_t [maxind1];
    r1[i] = new Float_t [maxind1];
    zc2[i] = new Float_t [maxind2];
    phi2[i] = new Float_t [maxind2];
    r2[i] = new Float_t [maxind2];
  }

  memset(ind1,0,nPhiBins*sizeof(Int_t));
  for(Int_t module= fFirstL1; module<=fLastL1;module++){
    if(module%4==0 || module%4==3)continue;
    branch->GetEvent(module);
    if(fUseV2Clusters){
      Clusters2RecPoints(clusters,module,itsRec);
    }
    Int_t nrecp1 = itsRec->GetEntries();
    for(Int_t j=0;j<nrecp1;j++){
      AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
      lc[0]=recp->GetX();
      lc[2]=recp->GetZ();
      geom->LtoG(module,lc,gc);
      gc[0]-=fX0;
      gc[1]-=fY0;
      Float_t xc1,yc1;
      xc1=gc[0];
      yc1=gc[1];
      Float_t phi = TMath::ATan2(gc[1],gc[0]);
      if(phi<0)phi=2*TMath::Pi()+phi;
      Int_t bin = (Int_t)(phi/phiBinSize);
      Int_t ind = ind1[bin];
      if(ind<maxind1) {
	phi1[bin][ind] = phi;
	zc1[bin][ind]=gc[2]/fStepFine;
	r1[bin][ind]=sqrt(xc1*xc1+yc1*yc1);
	ind1[bin]++;
      }
    }
    fITS->ResetRecPoints();
  }
  memset(ind2,0,nPhiBins*sizeof(Int_t));
  for(Int_t module= fFirstL2; module<=fLastL2;module++){
    branch->GetEvent(module);
    if(fUseV2Clusters){
      Clusters2RecPoints(clusters,module,itsRec);
    }
    Int_t nrecp2 = itsRec->GetEntries();
    for(Int_t j=0;j<nrecp2;j++){
      AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
      lc[0]=recp->GetX();
      lc[2]=recp->GetZ();
      geom->LtoG(module,lc,gc);
      gc[0]-=fX0;
      gc[1]-=fY0;
      Float_t xc2,yc2;
      xc2=gc[0];
      yc2=gc[1];
      Float_t phi = TMath::ATan2(gc[1],gc[0]);
      if(phi<0)phi=2*TMath::Pi()+phi;
      Int_t bin = (Int_t)(phi/phiBinSize+0.5);
      if(bin==nPhiBins) bin = 0;
      Int_t ind = ind2[bin];
      if(ind<maxind2) {
	phi2[bin][ind] = phi;
	zc2[bin][ind]=gc[2]/fStepFine;
	r2[bin][ind]=sqrt(xc2*xc2+yc2*yc2);
	ind2[bin]++;
      }
    }
    fITS->ResetRecPoints();
  }
  Int_t nbinfine = static_cast<Int_t>((fHighLim-fLowLim)/fStepFine);
  Float_t lowz = fLowLim/fStepFine;
  Float_t highz = fHighLim/fStepFine;
  Int_t *harray = new Int_t[nbinfine];
  memset(harray,0,nbinfine*sizeof(Int_t));
  for(Int_t ibin=0;ibin<nPhiBins;ibin++) {
    Float_t *pphi1 = phi1[ibin];
    Float_t *pr1 = r1[ibin];
    Float_t *pzc1 = zc1[ibin];
    for(Int_t i=0;i<ind1[ibin];i++){
      for(Int_t jbin=ibin;jbin<=(ibin+1);jbin++) {
	Int_t jbin2 = jbin;
	if(jbin==nPhiBins) jbin2 = 0;
	Float_t *pphi2 = phi2[jbin2];
	Float_t *pr2 = r2[jbin2];
	Float_t *pzc2 = zc2[jbin2];
	for(Int_t j=0;j<ind2[jbin2];j++){
	  Float_t diff = TMath::Abs(pphi2[j]-pphi1[i]);
	  if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
	  if(diff<fDiffPhiMax){
	    Float_t zr0=(pr2[j]*pzc1[i]-pr1[i]*pzc2[j])/(pr2[j]-pr1[i]);
	    if(zr0>lowz && zr0<highz) {
	      Int_t bin = (Int_t)(zr0-lowz);
	      harray[bin]++;
	    }
	  }
	}
      }
    }
  }
  for(Int_t i=0;i<nPhiBins;i++) {
    delete [] zc1[i];
    delete [] phi1[i];
    delete [] r1[i];
    delete [] zc2[i];
    delete [] phi2[i];
    delete [] r2[i];
  }
  delete [] zc1;
  delete [] phi1;
  delete [] r1;
  delete [] ind1;
  delete [] zc2;
  delete [] phi2;
  delete [] r2;
  delete [] ind2;

  Int_t nbinratio = (Int_t)(fStepCoarse/fStepFine+0.5);
  Int_t nbincoarse = nbinfine/nbinratio;

  if(fZCombc)delete fZCombc;
  fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);
  if(fZCombf)delete fZCombf;
  fZCombf = new TH1F("fZCombf","Z",nbinfine,fLowLim,fLowLim+nbinfine*fStepFine);
  Stat_t contents=0;
  Int_t counter=0;
  for(Int_t bin=0;bin<nbinfine;bin++) {
    fZCombf->SetBinContent(bin+1,(Stat_t)harray[bin]);
    fZCombf->SetBinError(bin+1,TMath::Sqrt((Stat_t)harray[bin]));
    contents+=(Stat_t)harray[bin];
    counter++;
    if(counter==nbinratio) {
      Int_t binc=bin/nbinratio; 
      fZCombc->SetBinContent(binc+1,contents);
      fZCombc->SetBinError(binc+1,TMath::Sqrt(contents));
      contents=0;
      counter=0;
    }
  }

  delete [] harray;

  if(fZCombc->GetEntries()==0){
    Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
    ResetHistograms();
    return fCurrentVertex;
  }
  //  else {
  //    cout<<"Number of entries in hist. "<<fZCombc->GetEntries()<<endl;
  //  }
  Int_t bi = fZCombc->GetMaximumBin();
  Float_t centre = fZCombc->GetBinCenter(bi);
  Int_t n1 = static_cast<Int_t>((centre-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
  Int_t n2 = static_cast<Int_t>((centre+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
  Int_t niter = 0;
  Bool_t goon = kTRUE;
  Int_t num;
  while(goon){
    fZFound = 0.;
    fZsig = 0.;
    num=0;
    for(Int_t n=n1;n<=n2;n++){
      fZFound+=fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
      num+=static_cast<Int_t>(fZCombf->GetBinContent(n));
      fZsig+=fZCombf->GetBinCenter(n)*fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
    }
    if(num<2){
      fZsig = 0.;
    }
    else {
      Float_t radi =  fZsig/(num-1)-fZFound*fZFound/num/(num-1);
      if(radi>0.)fZsig=TMath::Sqrt(radi);
      else fZsig=0.;
      fZFound/=num;
    }
    goon = TMath::Abs(TMath::Abs(fZFound-fZCombf->GetBinCenter(n1))-TMath::Abs(fZFound-fZCombf->GetBinCenter(n2)))>fTolerance;
    n1 = static_cast<Int_t>((fZFound-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
    n2 = static_cast<Int_t>((fZFound+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
    niter++;
    if(niter>=10){
      goon = kFALSE;
      Warning("FindVertexForCurrentEvent","The procedure dows not converge\n");
    }
  }
  //  cout<<"Numer of Iterations "<<niter<<endl<<endl;
  fCurrentVertex = new AliESDVertex(fZFound,fZsig,num);
  fCurrentVertex->SetTitle("vertexer: HLT");
  ResetHistograms();
  PrintStatus();
  return fCurrentVertex;
}
Commit	Line	Data
dd36288a	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15	#include "AliL3ITSVertexerZ.h"
	16	#include<TString.h>
	17	#include<TH1.h>
	18	#include<TMath.h>
	19	#include <AliRun.h>
	20	#include <AliITS.h>
	21	#include "AliITSLoader.h"
	22	#include <AliITSgeom.h>
	23	#include <AliITSRecPoint.h>
	24
	25	//-------------------------------------------------------------------------
	26	// Implementation of the HLT ITS vertexer class
	27	//
	28	// Origin: Cvetan Cheshkov, CERN, Cvetan.Cheshkov@cern.ch
	29	//-------------------------------------------------------------------------
	30
	31	ClassImp(AliL3ITSVertexerZ)
	32
	33	AliL3ITSVertexerZ::AliL3ITSVertexerZ(TString filename,Float_t x0, Float_t y0):AliITSVertexerZ(filename,x0,y0)
	34	{
	35	// Standard Constructor
	36	}
	37
	38	//______________________________________________________________________
	39	AliESDVertex* AliL3ITSVertexerZ::FindVertexForCurrentEvent(Int_t evnumber){
	40	// Defines the AliESDVertex for the current event
	41
	42	fCurrentVertex = 0;
	43	AliRunLoader *rl =AliRunLoader::GetRunLoader();
	44	AliITSLoader* itsLoader = (AliITSLoader*) rl->GetLoader("ITSLoader");
	45	itsLoader->LoadRecPoints();
	46	rl->GetEvent(evnumber);
	47
	48	if(!fITS) {
	49	fITS = (AliITS*)gAlice->GetModule("ITS");
	50	if(!fITS) {
	51	Error("FindVertexForCurrentEvent","AliITS object was not found");
	52	return fCurrentVertex;
	53	}
	54	}
	55
	56	fITS->SetTreeAddress();
	57
	58	AliITSgeom *geom = fITS->GetITSgeom();
	59
	60	TTree *tR = itsLoader->TreeR();
	61	TClonesArray *itsRec = 0;
	62	Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
	63	Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
	64	Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
65	Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
66
67	itsRec = fITS->RecPoints();
68
69	//cout<<"Address of itsRec = "<<itsRec<<endl;
70	TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
71	TBranch *branch;
72	if(fUseV2Clusters){
73	branch = tR->GetBranch("Clusters");
74	branch->SetAddress(&clusters);
75	}
76	else {
77	branch = tR->GetBranch("ITSRecPoints");
78	}
79
80
81	Int_t nrpL1 = 0;
82	Int_t nrpL2 = 0;
83	for(Int_t module= fFirstL1; module<=fLastL1;module++){
84	if(module%4==0 \|\| module%4==3)continue;
85	// cout<<"Procesing module "<<module<<" ";
86	branch->GetEvent(module);
87	// cout<<"Number of clusters "<<clusters->GetEntries()<<endl;
88	if(fUseV2Clusters){
89	Clusters2RecPoints(clusters,module,itsRec);
90	}
91	nrpL1+= itsRec->GetEntries();
92	fITS->ResetRecPoints();
93	}
94	for(Int_t module= fFirstL2; module<=fLastL2;module++){
95	branch->GetEvent(module);
96	if(fUseV2Clusters){
97	Clusters2RecPoints(clusters,module,itsRec);
98	}
99	nrpL2+= itsRec->GetEntries();
100	fITS->ResetRecPoints();
101	}
102	if(nrpL1 == 0 \|\| nrpL2 == 0){
103	ResetHistograms();
104	return fCurrentVertex;
105	}
106
107	Int_t nPhiBins = (Int_t)(TMath::Pi()/fDiffPhiMax);
108	Float_t phiBinSize = 2*TMath::Pi()/(Float_t)nPhiBins;
109
110	Int_t maxind1 = 2*nrpL1/nPhiBins;
111	Float_t *zc1 = new Float_t [nPhiBins];
112	Float_t *phi1 = new Float_t [nPhiBins];
113	Float_t *r1 = new Float_t [nPhiBins];
114	Int_t *ind1 = new Int_t [nPhiBins];
115	Int_t maxind2 = 2*nrpL2/nPhiBins;
116	Float_t *zc2 = new Float_t [nPhiBins];
117	Float_t *phi2 = new Float_t [nPhiBins];
118	Float_t *r2 = new Float_t [nPhiBins];
119	Int_t *ind2 = new Int_t [nPhiBins];
120	for(Int_t i=0;i<nPhiBins;i++) {
121	zc1[i] = new Float_t [maxind1];
122	phi1[i] = new Float_t [maxind1];
123	r1[i] = new Float_t [maxind1];
124	zc2[i] = new Float_t [maxind2];
125	phi2[i] = new Float_t [maxind2];
126	r2[i] = new Float_t [maxind2];
127	}
128
129	memset(ind1,0,nPhiBins*sizeof(Int_t));
130	for(Int_t module= fFirstL1; module<=fLastL1;module++){
131	if(module%4==0 \|\| module%4==3)continue;
132	branch->GetEvent(module);
133	if(fUseV2Clusters){
134	Clusters2RecPoints(clusters,module,itsRec);
135	}
136	Int_t nrecp1 = itsRec->GetEntries();
137	for(Int_t j=0;j<nrecp1;j++){
138	AliITSRecPoint recp = (AliITSRecPoint)itsRec->At(j);
139	lc[0]=recp->GetX();
140	lc[2]=recp->GetZ();
141	geom->LtoG(module,lc,gc);
142	gc[0]-=fX0;
143	gc[1]-=fY0;
144	Float_t xc1,yc1;
145	xc1=gc[0];
146	yc1=gc[1];
147	Float_t phi = TMath::ATan2(gc[1],gc[0]);
148	if(phi<0)phi=2*TMath::Pi()+phi;
149	Int_t bin = (Int_t)(phi/phiBinSize);
150	Int_t ind = ind1[bin];
151	if(ind<maxind1) {
152	phi1[bin][ind] = phi;
153	zc1[bin][ind]=gc[2]/fStepFine;
154	r1[bin][ind]=sqrt(xc1xc1+yc1yc1);
155	ind1[bin]++;
156	}
157	}
158	fITS->ResetRecPoints();
159	}
160	memset(ind2,0,nPhiBins*sizeof(Int_t));
161	for(Int_t module= fFirstL2; module<=fLastL2;module++){
162	branch->GetEvent(module);
163	if(fUseV2Clusters){
164	Clusters2RecPoints(clusters,module,itsRec);
165	}
166	Int_t nrecp2 = itsRec->GetEntries();
167	for(Int_t j=0;j<nrecp2;j++){
168	AliITSRecPoint recp = (AliITSRecPoint)itsRec->At(j);
169	lc[0]=recp->GetX();
170	lc[2]=recp->GetZ();
171	geom->LtoG(module,lc,gc);
172	gc[0]-=fX0;
173	gc[1]-=fY0;
174	Float_t xc2,yc2;
175	xc2=gc[0];
176	yc2=gc[1];
177	Float_t phi = TMath::ATan2(gc[1],gc[0]);
178	if(phi<0)phi=2*TMath::Pi()+phi;
179	Int_t bin = (Int_t)(phi/phiBinSize+0.5);
180	if(bin==nPhiBins) bin = 0;
181	Int_t ind = ind2[bin];
182	if(ind<maxind2) {
183	phi2[bin][ind] = phi;
184	zc2[bin][ind]=gc[2]/fStepFine;
185	r2[bin][ind]=sqrt(xc2xc2+yc2yc2);
186	ind2[bin]++;
187	}
188	}
189	fITS->ResetRecPoints();
190	}
191	Int_t nbinfine = static_cast<Int_t>((fHighLim-fLowLim)/fStepFine);
192	Float_t lowz = fLowLim/fStepFine;
193	Float_t highz = fHighLim/fStepFine;
194	Int_t *harray = new Int_t[nbinfine];
195	memset(harray,0,nbinfine*sizeof(Int_t));
196	for(Int_t ibin=0;ibin<nPhiBins;ibin++) {
197	Float_t *pphi1 = phi1[ibin];
198	Float_t *pr1 = r1[ibin];
199	Float_t *pzc1 = zc1[ibin];
200	for(Int_t i=0;i<ind1[ibin];i++){
201	for(Int_t jbin=ibin;jbin<=(ibin+1);jbin++) {
202	Int_t jbin2 = jbin;
203	if(jbin==nPhiBins) jbin2 = 0;
204	Float_t *pphi2 = phi2[jbin2];
205	Float_t *pr2 = r2[jbin2];
206	Float_t *pzc2 = zc2[jbin2];
207	for(Int_t j=0;j<ind2[jbin2];j++){
208	Float_t diff = TMath::Abs(pphi2[j]-pphi1[i]);
209	if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
210	if(diff<fDiffPhiMax){
211	Float_t zr0=(pr2[j]pzc1[i]-pr1[i]pzc2[j])/(pr2[j]-pr1[i]);
212	if(zr0>lowz && zr0<highz) {
213	Int_t bin = (Int_t)(zr0-lowz);
214	harray[bin]++;
215	}
216	}
217	}
218	}
219	}
220	}
221	for(Int_t i=0;i<nPhiBins;i++) {
222	delete [] zc1[i];
223	delete [] phi1[i];
224	delete [] r1[i];
225	delete [] zc2[i];
226	delete [] phi2[i];
227	delete [] r2[i];
228	}
229	delete [] zc1;
230	delete [] phi1;
231	delete [] r1;
232	delete [] ind1;
233	delete [] zc2;
234	delete [] phi2;
235	delete [] r2;
236	delete [] ind2;
237
238	Int_t nbinratio = (Int_t)(fStepCoarse/fStepFine+0.5);
239	Int_t nbincoarse = nbinfine/nbinratio;
240
241	if(fZCombc)delete fZCombc;
242	fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);
243	if(fZCombf)delete fZCombf;
244	fZCombf = new TH1F("fZCombf","Z",nbinfine,fLowLim,fLowLim+nbinfine*fStepFine);
245	Stat_t contents=0;
246	Int_t counter=0;
247	for(Int_t bin=0;bin<nbinfine;bin++) {
248	fZCombf->SetBinContent(bin+1,(Stat_t)harray[bin]);
249	fZCombf->SetBinError(bin+1,TMath::Sqrt((Stat_t)harray[bin]));
250	contents+=(Stat_t)harray[bin];
251	counter++;
252	if(counter==nbinratio) {
253	Int_t binc=bin/nbinratio;
254	fZCombc->SetBinContent(binc+1,contents);
255	fZCombc->SetBinError(binc+1,TMath::Sqrt(contents));
256	contents=0;
257	counter=0;
258	}
259	}
260
261	delete [] harray;
262
263	if(fZCombc->GetEntries()==0){
264	Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
265	ResetHistograms();
266	return fCurrentVertex;
267	}
268	// else {
269	// cout<<"Number of entries in hist. "<<fZCombc->GetEntries()<<endl;
270	// }
271	Int_t bi = fZCombc->GetMaximumBin();
272	Float_t centre = fZCombc->GetBinCenter(bi);
273	Int_t n1 = static_cast<Int_t>((centre-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
274	Int_t n2 = static_cast<Int_t>((centre+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
275	Int_t niter = 0;
276	Bool_t goon = kTRUE;
277	Int_t num;
278	while(goon){
279	fZFound = 0.;
280	fZsig = 0.;
281	num=0;
282	for(Int_t n=n1;n<=n2;n++){
283	fZFound+=fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
284	num+=static_cast<Int_t>(fZCombf->GetBinContent(n));
285	fZsig+=fZCombf->GetBinCenter(n)fZCombf->GetBinCenter(n)fZCombf->GetBinContent(n);
286	}
287	if(num<2){
288	fZsig = 0.;
289	}
290	else {
291	Float_t radi = fZsig/(num-1)-fZFound*fZFound/num/(num-1);
292	if(radi>0.)fZsig=TMath::Sqrt(radi);
293	else fZsig=0.;
294	fZFound/=num;
295	}
296	goon = TMath::Abs(TMath::Abs(fZFound-fZCombf->GetBinCenter(n1))-TMath::Abs(fZFound-fZCombf->GetBinCenter(n2)))>fTolerance;
297	n1 = static_cast<Int_t>((fZFound-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
298	n2 = static_cast<Int_t>((fZFound+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
299	niter++;
300	if(niter>=10){
301	goon = kFALSE;
302	Warning("FindVertexForCurrentEvent","The procedure dows not converge\n");
303	}
304	}
305	// cout<<"Numer of Iterations "<<niter<<endl<<endl;
306	fCurrentVertex = new AliESDVertex(fZFound,fZsig,num);
307	fCurrentVertex->SetTitle("vertexer: HLT");
308	ResetHistograms();
309	PrintStatus();
310	return fCurrentVertex;
311	}