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

// $Id$

/**************************************************************************
 * 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 "AliHLTITSVertexerZ.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>
#include <AliITSclusterV2.h>

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

ClassImp(AliHLTITSVertexerZ)

AliHLTITSVertexerZ::AliHLTITSVertexerZ():
  AliITSVertexerZ(),
  fZCombf(0),
  fStepFine(0)
{
  // Constructor in case that there is no runloader
  SetBinWidthFine();
}

AliHLTITSVertexerZ::AliHLTITSVertexerZ(TString filename,Float_t x0, Float_t y0):
  AliITSVertexerZ(filename,x0,y0),
  fZCombf(0),
  fStepFine(0)
{
  // Standard Constructor
  SetBinWidthFine();
}

AliHLTITSVertexerZ::~AliHLTITSVertexerZ()
{
  // Destructor
  if (fZCombf) delete fZCombf;
}

//______________________________________________________________________
AliHLTITSVertexerZ::AliHLTITSVertexerZ(const AliHLTITSVertexerZ &vtxr) :
  AliITSVertexerZ(vtxr),
  fZCombf(vtxr.fZCombf),
  fStepFine(vtxr.fStepFine)
{
  // Copy constructor

}

//______________________________________________________________________
AliHLTITSVertexerZ& AliHLTITSVertexerZ::operator=(const AliHLTITSVertexerZ&  vtxr )
{
  // Assignment operator
  this->~AliHLTITSVertexerZ();
  new(this) AliHLTITSVertexerZ(vtxr);
  fZCombf = vtxr.fZCombf;
  fStepFine = vtxr.fStepFine;

  return *this;
}
//______________________________________________________________________
AliESDVertex* AliHLTITSVertexerZ::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();
  rl->CdGAFile();

  AliITSgeom* geom=itsLoader->GetITSgeom();
  TTree *tR = itsLoader->TreeR();
  
  return FindVertexForCurrentEvent(geom,tR);
}

//______________________________________________________________________
AliESDVertex* AliHLTITSVertexerZ::FindVertexForCurrentEvent(AliITSgeom *geom,TTree *tR){
  // Defines the AliESDVertex for the current event

  fCurrentVertex = 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.;

  TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
  TBranch *branch;
  branch = tR->GetBranch("Clusters");
  branch->SetAddress(&clusters);

  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<<" ";
    tR->GetEvent(module);
    //    cout<<"Number of clusters "<<clusters->GetEntries()<<endl;
    nrpL1+= clusters->GetEntriesFast();
  }
  for(Int_t module= fFirstL2; module<=fLastL2;module++){
    tR->GetEvent(module);
    nrpL2+= clusters->GetEntriesFast();
  }
  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];
  }
  
  Float_t yshift = 0;
  Float_t zshift[4] = {-10.708000, -3.536000, 3.536000, 10.708000};

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