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

/**************************************************************************
 * Copyright(c) 1998-2003, 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 <AliITSVertexerZ.h>
#include <TString.h>
#include "AliITSLoader.h"
#include<TBranch.h>
#include<TClonesArray.h>
#include<TH1.h>
#include<TTree.h>
#include <AliITSgeom.h>
#include "AliITSDetTypeRec.h"
#include <AliITSRecPoint.h>

/////////////////////////////////////////////////////////////////
// this class implements a fast method to determine
// the Z coordinate of the primary vertex
// for p-p collisions it seems to give comparable or better results
// with respect to what obtained with AliITSVertexerPPZ
// It can be used successfully with Pb-Pb collisions
////////////////////////////////////////////////////////////////

ClassImp(AliITSVertexerZ)


//______________________________________________________________________
AliITSVertexerZ::AliITSVertexerZ():AliITSVertexer() {
  // Default constructor
  SetDiffPhiMax(0);
  fX0 = 0.;
  fY0 = 0.;
  SetFirstLayerModules(0);
  SetSecondLayerModules(0);
  fZFound = 0;
  fZsig = 0.;
  fZCombc = 0;
  fZCombf = 0;
  SetLowLimit(0.);
  SetHighLimit(0.);
  SetBinWidthCoarse(0.);
  SetBinWidthFine(0.);
  SetTolerance(0.);
}

//______________________________________________________________________
AliITSVertexerZ::AliITSVertexerZ(TString fn, Float_t x0, Float_t y0):AliITSVertexer(fn) {
  // Standard Constructor
  SetDiffPhiMax();
  fX0 = x0;
  fY0 = y0;
  SetFirstLayerModules();
  SetSecondLayerModules();
  fZFound = 0;
  fZsig = 0.;
  fZCombc = 0;
  fZCombf = 0;
  SetLowLimit();
  SetHighLimit();
  SetBinWidthCoarse();
  SetBinWidthFine();
  SetTolerance();

}

//______________________________________________________________________
AliITSVertexerZ::AliITSVertexerZ(const AliITSVertexerZ &vtxr) : AliITSVertexer(vtxr) {
  // Copy constructor
  // Copies are not allowed. The method is protected to avoid misuse.
  Error("AliITSVertexerZ","Copy constructor not allowed\n");
}

//______________________________________________________________________
AliITSVertexerZ& AliITSVertexerZ::operator=(const AliITSVertexerZ& /* vtxr */){
  // Assignment operator
  // Assignment is not allowed. The method is protected to avoid misuse.
  Error("= operator","Assignment operator not allowed\n");
  return *this;
}


//______________________________________________________________________
AliITSVertexerZ::~AliITSVertexerZ() {
  // Default Destructor
  //fITS = 0;
  if(fZCombc)delete fZCombc;
  if(fZCombf)delete fZCombf;
}

//______________________________________________________________________
AliESDVertex* AliITSVertexerZ::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* geom = (AliITSgeom*)gDirectory->Get("AliITSgeom");
  olddir->cd();

  itsLoader->LoadRecPoints();
  rl->GetEvent(evnumber);

  AliITSDetTypeRec detTypeRec;

  TTree *tR = itsLoader->TreeR();
  detTypeRec.SetTreeAddressR(tR);
  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 = detTypeRec.RecPoints();
  TBranch *branch;
  branch = tR->GetBranch("ITSRecPoints");

  Int_t nbinfine = static_cast<Int_t>((fHighLim-fLowLim)/fStepFine);
  Int_t nbincoarse = static_cast<Int_t>((fHighLim-fLowLim)/fStepCoarse);
  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);

  Int_t nrpL1 = 0;
  Int_t nrpL2 = 0;
  for(Int_t module= fFirstL1; module<=fLastL1;module++){
    if(module%4==0 || module%4==3)continue;
    branch->GetEvent(module);
    nrpL1+= itsRec->GetEntries();
    detTypeRec.ResetRecPoints();
  }
  for(Int_t module= fFirstL2; module<=fLastL2;module++){
    branch->GetEvent(module);
    nrpL2+= itsRec->GetEntries();
    detTypeRec.ResetRecPoints();
  }
  if(nrpL1 == 0 || nrpL2 == 0){
    ResetHistograms();
    return fCurrentVertex;
  }
  Float_t *xc1 = new Float_t [nrpL1];
  Float_t *yc1 = new Float_t [nrpL1];
  Float_t *zc1 = new Float_t [nrpL1];
  Float_t *phi1 = new Float_t [nrpL1];
  Float_t *xc2 = new Float_t [nrpL2];
  Float_t *yc2 = new Float_t [nrpL2];
  Float_t *zc2 = new Float_t [nrpL2];
  Float_t *phi2 = new Float_t [nrpL2];
  Int_t ind = 0;
  for(Int_t module= fFirstL1; module<=fLastL1;module++){
    if(module%4==0 || module%4==3)continue;
    branch->GetEvent(module);
    Int_t nrecp1 = itsRec->GetEntries();
    for(Int_t j=0;j<nrecp1;j++){
      AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
      lc[0]=recp->GetDetLocalX();
      lc[2]=recp->GetDetLocalZ();
      geom->LtoG(module,lc,gc);
      gc[0]-=fX0;
      gc[1]-=fY0;
      xc1[ind]=gc[0];
      yc1[ind]=gc[1];
      zc1[ind]=gc[2];
      phi1[ind] = TMath::ATan2(gc[1],gc[0]);
      if(phi1[ind]<0)phi1[ind]=2*TMath::Pi()+phi1[ind];
      ind++;
    }
    detTypeRec.ResetRecPoints();
  }
  ind = 0;
  for(Int_t module= fFirstL2; module<=fLastL2;module++){
    branch->GetEvent(module);
    Int_t nrecp2 = itsRec->GetEntries();
    for(Int_t j=0;j<nrecp2;j++){
      AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
      lc[0]=recp->GetDetLocalX();
      lc[2]=recp->GetDetLocalZ();
      geom->LtoG(module,lc,gc);
      gc[0]-=fX0;
      gc[1]-=fY0;
      xc2[ind]=gc[0];
      yc2[ind]=gc[1];
      zc2[ind]=gc[2];
      phi2[ind] = TMath::ATan2(gc[1],gc[0]);
      if(phi2[ind]<0)phi2[ind]=2*TMath::Pi()+phi2[ind];
      ind++;
    }
    detTypeRec.ResetRecPoints();
  }
  for(Int_t i=0;i<nrpL1;i++){
    Float_t r1=TMath::Sqrt(xc1[i]*xc1[i]+yc1[i]*yc1[i]);
    for(Int_t j=0;j<nrpL2;j++){
      Float_t diff = TMath::Abs(phi2[j]-phi1[i]);
      if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
      if(diff<fDiffPhiMax){
	Float_t r2=TMath::Sqrt(xc2[j]*xc2[j]+yc2[j]*yc2[j]);
	Float_t zr0=(r2*zc1[i]-r1*zc2[j])/(r2-r1);
	fZCombf->Fill(zr0);
	fZCombc->Fill(zr0);
      }
    }
  }
  delete [] xc1;
  delete [] yc1;
  delete [] zc1;
  delete [] phi1;
  delete [] xc2;
  delete [] yc2;
  delete [] zc2;
  delete [] phi2;
  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: B");
  ResetHistograms();
  return fCurrentVertex;
}

//_____________________________________________________________________
void AliITSVertexerZ::ResetHistograms(){
  // delete TH1 data members
  if(fZCombc)delete fZCombc;
  if(fZCombf)delete fZCombf;
  fZCombc = 0;
  fZCombf = 0;
}

//______________________________________________________________________
void AliITSVertexerZ::FindVertices(){
  // computes the vertices of the events in the range FirstEvent - LastEvent
  AliRunLoader *rl = AliRunLoader::GetRunLoader();
  AliITSLoader* itsLoader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
  itsLoader->ReloadRecPoints();
  for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
    cout<<"Processing event "<<i<<endl;
    rl->GetEvent(i);
    FindVertexForCurrentEvent(i);
    if(fCurrentVertex){
      WriteCurrentVertex();
    }
    else {
      if(fDebug>0){
	cout<<"Vertex not found for event "<<i<<endl;
	cout<<"fZFound = "<<fZFound<<", fZsig= "<<fZsig<<endl;
      }
    }
  }
}

//________________________________________________________
void AliITSVertexerZ::PrintStatus() const {
  // Print current status
  cout <<"=======================================================\n";
  cout <<" First layer first and last modules: "<<fFirstL1<<", ";
  cout <<fLastL1<<endl;
  cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
  cout <<fLastL2<<endl;
  cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
  cout <<"Limits for Z histograms: "<<fLowLim<<"; "<<fHighLim<<endl;
  cout <<"Bin sizes for coarse and fine z histos "<<fStepCoarse<<"; "<<fStepFine<<endl;
  cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
  cout <<" Debug flag: "<<fDebug<<endl;
  cout <<"First event to be processed "<<fFirstEvent;
  cout <<"\n Last event to be processed "<<fLastEvent<<endl;
 
 
  cout <<"=======================================================\n";
}
Commit	Line	Data
0c6af5c9	1	/**************************************************************************
	2	* Copyright(c) 1998-2003, 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 <AliITSVertexerZ.h>
0c6af5c9	16	#include <TString.h>
0c6af5c9	17	#include "AliITSLoader.h"
0c6af5c9	18	#include<TBranch.h>
	19	#include<TClonesArray.h>
	20	#include<TH1.h>
0c6af5c9	21	#include<TTree.h>
0c6af5c9	22	#include <AliITSgeom.h>
7d62fb64	23	#include "AliITSDetTypeRec.h"
0c6af5c9	24	#include <AliITSRecPoint.h>
	25
	26	/////////////////////////////////////////////////////////////////
	27	// this class implements a fast method to determine
	28	// the Z coordinate of the primary vertex
	29	// for p-p collisions it seems to give comparable or better results
	30	// with respect to what obtained with AliITSVertexerPPZ
	31	// It can be used successfully with Pb-Pb collisions
	32	////////////////////////////////////////////////////////////////
	33
	34	ClassImp(AliITSVertexerZ)
	35
	36
	37
	38	//______________________________________________________________________
	39	AliITSVertexerZ::AliITSVertexerZ():AliITSVertexer() {
	40	// Default constructor
	41	SetDiffPhiMax(0);
	42	fX0 = 0.;
	43	fY0 = 0.;
	44	SetFirstLayerModules(0);
	45	SetSecondLayerModules(0);
	46	fZFound = 0;
	47	fZsig = 0.;
0c6af5c9	48	fZCombc = 0;
	49	fZCombf = 0;
	50	SetLowLimit(0.);
	51	SetHighLimit(0.);
	52	SetBinWidthCoarse(0.);
	53	SetBinWidthFine(0.);
	54	SetTolerance(0.);
	55	}
	56
	57	//______________________________________________________________________
	58	AliITSVertexerZ::AliITSVertexerZ(TString fn, Float_t x0, Float_t y0):AliITSVertexer(fn) {
	59	// Standard Constructor
	60	SetDiffPhiMax();
	61	fX0 = x0;
	62	fY0 = y0;
	63	SetFirstLayerModules();
	64	SetSecondLayerModules();
	65	fZFound = 0;
	66	fZsig = 0.;
0c6af5c9	67	fZCombc = 0;
	68	fZCombf = 0;
	69	SetLowLimit();
	70	SetHighLimit();
	71	SetBinWidthCoarse();
	72	SetBinWidthFine();
	73	SetTolerance();
	74
	75	}
	76
	77	//______________________________________________________________________
	78	AliITSVertexerZ::AliITSVertexerZ(const AliITSVertexerZ &vtxr) : AliITSVertexer(vtxr) {
	79	// Copy constructor
	80	// Copies are not allowed. The method is protected to avoid misuse.
	81	Error("AliITSVertexerZ","Copy constructor not allowed\n");
	82	}
	83
	84	//______________________________________________________________________
	85	AliITSVertexerZ& AliITSVertexerZ::operator=(const AliITSVertexerZ& /* vtxr */){
	86	// Assignment operator
	87	// Assignment is not allowed. The method is protected to avoid misuse.
	88	Error("= operator","Assignment operator not allowed\n");
	89	return *this;
	90	}
	91
	92
	93	//______________________________________________________________________
	94	AliITSVertexerZ::~AliITSVertexerZ() {
	95	// Default Destructor
7d62fb64	96	//fITS = 0;
0c6af5c9	97	if(fZCombc)delete fZCombc;
	98	if(fZCombf)delete fZCombf;
	99	}
	100
	101	//______________________________________________________________________
d681bb2d	102	AliESDVertex* AliITSVertexerZ::FindVertexForCurrentEvent(Int_t evnumber){
d681bb2d	103	// Defines the AliESDVertex for the current event
7d62fb64	104
0c6af5c9	105	fCurrentVertex = 0;
0c6af5c9	106	AliRunLoader *rl =AliRunLoader::GetRunLoader();
7d62fb64	107	AliITSLoader* itsLoader = (AliITSLoader*)rl->GetLoader("ITSLoader");
	108	TDirectory * olddir = gDirectory;
	109	rl->CdGAFile();
	110	AliITSgeom* geom = (AliITSgeom*)gDirectory->Get("AliITSgeom");
	111	olddir->cd();
	112
2257f27e	113	itsLoader->LoadRecPoints();
0c6af5c9	114	rl->GetEvent(evnumber);
0c6af5c9	115
7d62fb64	116	AliITSDetTypeRec detTypeRec;
0c6af5c9	117
0c6af5c9	118	TTree *tR = itsLoader->TreeR();
7d62fb64	119	detTypeRec.SetTreeAddressR(tR);
0c6af5c9	120	TClonesArray *itsRec = 0;
	121	Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
	122	Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
	123	Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
	124	Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
	125
7d62fb64	126	itsRec = detTypeRec.RecPoints();
2257f27e	127	TBranch *branch;
00a7cc50	128	branch = tR->GetBranch("ITSRecPoints");
0c6af5c9	129
	130	Int_t nbinfine = static_cast<Int_t>((fHighLim-fLowLim)/fStepFine);
	131	Int_t nbincoarse = static_cast<Int_t>((fHighLim-fLowLim)/fStepCoarse);
	132	if(fZCombc)delete fZCombc;
	133	fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);
	134	if(fZCombf)delete fZCombf;
	135	fZCombf = new TH1F("fZCombf","Z",nbinfine,fLowLim,fLowLim+nbinfine*fStepFine);
	136
	137	Int_t nrpL1 = 0;
	138	Int_t nrpL2 = 0;
	139	for(Int_t module= fFirstL1; module<=fLastL1;module++){
	140	if(module%4==0 \|\| module%4==3)continue;
	141	branch->GetEvent(module);
	142	nrpL1+= itsRec->GetEntries();
7d62fb64	143	detTypeRec.ResetRecPoints();
0c6af5c9	144	}
	145	for(Int_t module= fFirstL2; module<=fLastL2;module++){
	146	branch->GetEvent(module);
	147	nrpL2+= itsRec->GetEntries();
7d62fb64	148	detTypeRec.ResetRecPoints();
0c6af5c9	149	}
0c6af5c9	150	if(nrpL1 == 0 \|\| nrpL2 == 0){
0c6af5c9	151	ResetHistograms();
	152	return fCurrentVertex;
	153	}
	154	Float_t *xc1 = new Float_t [nrpL1];
	155	Float_t *yc1 = new Float_t [nrpL1];
	156	Float_t *zc1 = new Float_t [nrpL1];
	157	Float_t *phi1 = new Float_t [nrpL1];
	158	Float_t *xc2 = new Float_t [nrpL2];
	159	Float_t *yc2 = new Float_t [nrpL2];
	160	Float_t *zc2 = new Float_t [nrpL2];
	161	Float_t *phi2 = new Float_t [nrpL2];
	162	Int_t ind = 0;
	163	for(Int_t module= fFirstL1; module<=fLastL1;module++){
	164	if(module%4==0 \|\| module%4==3)continue;
	165	branch->GetEvent(module);
	166	Int_t nrecp1 = itsRec->GetEntries();
	167	for(Int_t j=0;j<nrecp1;j++){
	168	AliITSRecPoint recp = (AliITSRecPoint)itsRec->At(j);
00a7cc50	169	lc[0]=recp->GetDetLocalX();
00a7cc50	170	lc[2]=recp->GetDetLocalZ();
0c6af5c9	171	geom->LtoG(module,lc,gc);
	172	gc[0]-=fX0;
	173	gc[1]-=fY0;
	174	xc1[ind]=gc[0];
	175	yc1[ind]=gc[1];
	176	zc1[ind]=gc[2];
	177	phi1[ind] = TMath::ATan2(gc[1],gc[0]);
	178	if(phi1[ind]<0)phi1[ind]=2*TMath::Pi()+phi1[ind];
	179	ind++;
	180	}
7d62fb64	181	detTypeRec.ResetRecPoints();
0c6af5c9	182	}
	183	ind = 0;
	184	for(Int_t module= fFirstL2; module<=fLastL2;module++){
	185	branch->GetEvent(module);
	186	Int_t nrecp2 = itsRec->GetEntries();
	187	for(Int_t j=0;j<nrecp2;j++){
	188	AliITSRecPoint recp = (AliITSRecPoint)itsRec->At(j);
00a7cc50	189	lc[0]=recp->GetDetLocalX();
00a7cc50	190	lc[2]=recp->GetDetLocalZ();
0c6af5c9	191	geom->LtoG(module,lc,gc);
	192	gc[0]-=fX0;
	193	gc[1]-=fY0;
	194	xc2[ind]=gc[0];
	195	yc2[ind]=gc[1];
	196	zc2[ind]=gc[2];
	197	phi2[ind] = TMath::ATan2(gc[1],gc[0]);
	198	if(phi2[ind]<0)phi2[ind]=2*TMath::Pi()+phi2[ind];
	199	ind++;
	200	}
7d62fb64	201	detTypeRec.ResetRecPoints();
0c6af5c9	202	}
	203	for(Int_t i=0;i<nrpL1;i++){
	204	Float_t r1=TMath::Sqrt(xc1[i]xc1[i]+yc1[i]yc1[i]);
	205	for(Int_t j=0;j<nrpL2;j++){
	206	Float_t diff = TMath::Abs(phi2[j]-phi1[i]);
	207	if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
	208	if(diff<fDiffPhiMax){
	209	Float_t r2=TMath::Sqrt(xc2[j]xc2[j]+yc2[j]yc2[j]);
	210	Float_t zr0=(r2zc1[i]-r1zc2[j])/(r2-r1);
	211	fZCombf->Fill(zr0);
	212	fZCombc->Fill(zr0);
	213	}
	214	}
	215	}
	216	delete [] xc1;
	217	delete [] yc1;
	218	delete [] zc1;
	219	delete [] phi1;
	220	delete [] xc2;
	221	delete [] yc2;
	222	delete [] zc2;
	223	delete [] phi2;
	224	if(fZCombc->GetEntries()==0){
	225	Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
	226	ResetHistograms();
	227	return fCurrentVertex;
	228	}
2257f27e	229	// else {
	230	// cout<<"Number of entries in hist. "<<fZCombc->GetEntries()<<endl;
	231	// }
0c6af5c9	232	Int_t bi = fZCombc->GetMaximumBin();
	233	Float_t centre = fZCombc->GetBinCenter(bi);
	234	Int_t n1 = static_cast<Int_t>((centre-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
	235	Int_t n2 = static_cast<Int_t>((centre+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
	236	Int_t niter = 0;
	237	Bool_t goon = kTRUE;
	238	Int_t num;
	239	while(goon){
	240	fZFound = 0.;
	241	fZsig = 0.;
	242	num=0;
	243	for(Int_t n=n1;n<=n2;n++){
	244	fZFound+=fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
	245	num+=static_cast<Int_t>(fZCombf->GetBinContent(n));
	246	fZsig+=fZCombf->GetBinCenter(n)fZCombf->GetBinCenter(n)fZCombf->GetBinContent(n);
	247	}
	248	if(num<2){
0c6af5c9	249	fZsig = 0.;
	250	}
	251	else {
979e3647	252	Float_t radi = fZsig/(num-1)-fZFound*fZFound/num/(num-1);
	253	if(radi>0.)fZsig=TMath::Sqrt(radi);
	254	else fZsig=0.;
0c6af5c9	255	fZFound/=num;
	256	}
	257	goon = TMath::Abs(TMath::Abs(fZFound-fZCombf->GetBinCenter(n1))-TMath::Abs(fZFound-fZCombf->GetBinCenter(n2)))>fTolerance;
0c6af5c9	258	n1 = static_cast<Int_t>((fZFound-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
	259	n2 = static_cast<Int_t>((fZFound+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
	260	niter++;
	261	if(niter>=10){
	262	goon = kFALSE;
	263	Warning("FindVertexForCurrentEvent","The procedure dows not converge\n");
	264	}
	265	}
2257f27e	266	// cout<<"Numer of Iterations "<<niter<<endl<<endl;
d681bb2d	267	fCurrentVertex = new AliESDVertex(fZFound,fZsig,num);
0c6af5c9	268	fCurrentVertex->SetTitle("vertexer: B");
	269	ResetHistograms();
	270	return fCurrentVertex;
	271	}
	272
	273	//_____________________________________________________________________
	274	void AliITSVertexerZ::ResetHistograms(){
	275	// delete TH1 data members
	276	if(fZCombc)delete fZCombc;
	277	if(fZCombf)delete fZCombf;
	278	fZCombc = 0;
	279	fZCombf = 0;
	280	}
	281
	282	//______________________________________________________________________
	283	void AliITSVertexerZ::FindVertices(){
	284	// computes the vertices of the events in the range FirstEvent - LastEvent
	285	AliRunLoader *rl = AliRunLoader::GetRunLoader();
	286	AliITSLoader* itsLoader = (AliITSLoader*) rl->GetLoader("ITSLoader");
	287	itsLoader->ReloadRecPoints();
	288	for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
	289	cout<<"Processing event "<<i<<endl;
	290	rl->GetEvent(i);
	291	FindVertexForCurrentEvent(i);
	292	if(fCurrentVertex){
	293	WriteCurrentVertex();
	294	}
	295	else {
	296	if(fDebug>0){
	297	cout<<"Vertex not found for event "<<i<<endl;
	298	cout<<"fZFound = "<<fZFound<<", fZsig= "<<fZsig<<endl;
	299	}
	300	}
	301	}
	302	}
	303
	304	//________________________________________________________
	305	void AliITSVertexerZ::PrintStatus() const {
	306	// Print current status
	307	cout <<"=======================================================\n";
	308	cout <<" First layer first and last modules: "<<fFirstL1<<", ";
	309	cout <<fLastL1<<endl;
	310	cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
	311	cout <<fLastL2<<endl;
	312	cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
	313	cout <<"Limits for Z histograms: "<<fLowLim<<"; "<<fHighLim<<endl;
	314	cout <<"Bin sizes for coarse and fine z histos "<<fStepCoarse<<"; "<<fStepFine<<endl;
	315	cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
	316	cout <<" Debug flag: "<<fDebug<<endl;
	317	cout <<"First event to be processed "<<fFirstEvent;
	318	cout <<"\n Last event to be processed "<<fLastEvent<<endl;
	319
	320
	321	cout <<"=======================================================\n";
	322	}
	323