[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 <Riostream.h>
#include <TString.h>
#include <AliITS.h>
#include "AliITSLoader.h"
#include <AliRun.h>
#include<TBranch.h>
#include<TClonesArray.h>
#include<TH1.h>
#include<TMath.h>
#include<TTree.h>
#include <AliITSgeom.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.;
  fITS = 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.;
  fITS = 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");
  itsLoader->ReloadRecPoints();

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