[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->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 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;
    //   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;
  }
  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);
    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;
      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);
    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;
      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;
    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");
2257f27e	113	itsLoader->LoadRecPoints();
0c6af5c9	114	rl->GetEvent(evnumber);
	115
	116	if(!fITS) {
	117	fITS = (AliITS*)gAlice->GetModule("ITS");
	118	if(!fITS) {
	119	Error("FindVertexForCurrentEvent","AliITS object was not found");
	120	return fCurrentVertex;
	121	}
	122	}
	123
	124	fITS->SetTreeAddress();
	125
	126	AliITSgeom *geom = fITS->GetITSgeom();
	127
	128	TTree *tR = itsLoader->TreeR();
	129	TClonesArray *itsRec = 0;
	130	Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
	131	Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
	132	Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
	133	Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
	134
2257f27e	135	itsRec = fITS->RecPoints();
	136
	137	//cout<<"Address of itsRec = "<<itsRec<<endl;
	138	TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
	139	TBranch *branch;
	140	if(fUseV2Clusters){
	141	branch = tR->GetBranch("Clusters");
	142	branch->SetAddress(&clusters);
	143	}
	144	else {
	145	branch = tR->GetBranch("ITSRecPoints");
	146	}
0c6af5c9	147
	148	Int_t nbinfine = static_cast<Int_t>((fHighLim-fLowLim)/fStepFine);
	149	Int_t nbincoarse = static_cast<Int_t>((fHighLim-fLowLim)/fStepCoarse);
	150	if(fZCombc)delete fZCombc;
	151	fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);
	152	if(fZCombf)delete fZCombf;
	153	fZCombf = new TH1F("fZCombf","Z",nbinfine,fLowLim,fLowLim+nbinfine*fStepFine);
	154
	155	Int_t nrpL1 = 0;
	156	Int_t nrpL2 = 0;
	157	for(Int_t module= fFirstL1; module<=fLastL1;module++){
	158	if(module%4==0 \|\| module%4==3)continue;
2257f27e	159	// cout<<"Procesing module "<<module<<" ";
0c6af5c9	160	branch->GetEvent(module);
2257f27e	161	// cout<<"Number of clusters "<<clusters->GetEntries()<<endl;
	162	if(fUseV2Clusters){
	163	Clusters2RecPoints(clusters,module,itsRec);
	164	}
0c6af5c9	165	nrpL1+= itsRec->GetEntries();
	166	fITS->ResetRecPoints();
	167	}
	168	for(Int_t module= fFirstL2; module<=fLastL2;module++){
	169	branch->GetEvent(module);
2257f27e	170	if(fUseV2Clusters){
	171	Clusters2RecPoints(clusters,module,itsRec);
	172	}
0c6af5c9	173	nrpL2+= itsRec->GetEntries();
	174	fITS->ResetRecPoints();
	175	}
0c6af5c9	176	if(nrpL1 == 0 \|\| nrpL2 == 0){
0c6af5c9	177	ResetHistograms();
	178	return fCurrentVertex;
	179	}
	180	Float_t *xc1 = new Float_t [nrpL1];
	181	Float_t *yc1 = new Float_t [nrpL1];
	182	Float_t *zc1 = new Float_t [nrpL1];
	183	Float_t *phi1 = new Float_t [nrpL1];
	184	Float_t *xc2 = new Float_t [nrpL2];
	185	Float_t *yc2 = new Float_t [nrpL2];
	186	Float_t *zc2 = new Float_t [nrpL2];
	187	Float_t *phi2 = new Float_t [nrpL2];
	188	Int_t ind = 0;
	189	for(Int_t module= fFirstL1; module<=fLastL1;module++){
	190	if(module%4==0 \|\| module%4==3)continue;
	191	branch->GetEvent(module);
2257f27e	192	if(fUseV2Clusters){
	193	Clusters2RecPoints(clusters,module,itsRec);
	194	}
0c6af5c9	195	Int_t nrecp1 = itsRec->GetEntries();
	196	for(Int_t j=0;j<nrecp1;j++){
	197	AliITSRecPoint recp = (AliITSRecPoint)itsRec->At(j);
	198	lc[0]=recp->GetX();
	199	lc[2]=recp->GetZ();
	200	geom->LtoG(module,lc,gc);
	201	gc[0]-=fX0;
	202	gc[1]-=fY0;
	203	xc1[ind]=gc[0];
	204	yc1[ind]=gc[1];
	205	zc1[ind]=gc[2];
	206	phi1[ind] = TMath::ATan2(gc[1],gc[0]);
	207	if(phi1[ind]<0)phi1[ind]=2*TMath::Pi()+phi1[ind];
	208	ind++;
	209	}
	210	fITS->ResetRecPoints();
	211	}
	212	ind = 0;
	213	for(Int_t module= fFirstL2; module<=fLastL2;module++){
	214	branch->GetEvent(module);
2257f27e	215	if(fUseV2Clusters){
	216	Clusters2RecPoints(clusters,module,itsRec);
	217	}
0c6af5c9	218	Int_t nrecp2 = itsRec->GetEntries();
	219	for(Int_t j=0;j<nrecp2;j++){
	220	AliITSRecPoint recp = (AliITSRecPoint)itsRec->At(j);
	221	lc[0]=recp->GetX();
	222	lc[2]=recp->GetZ();
	223	geom->LtoG(module,lc,gc);
	224	gc[0]-=fX0;
	225	gc[1]-=fY0;
	226	xc2[ind]=gc[0];
	227	yc2[ind]=gc[1];
	228	zc2[ind]=gc[2];
	229	phi2[ind] = TMath::ATan2(gc[1],gc[0]);
	230	if(phi2[ind]<0)phi2[ind]=2*TMath::Pi()+phi2[ind];
	231	ind++;
	232	}
	233	fITS->ResetRecPoints();
	234	}
	235	for(Int_t i=0;i<nrpL1;i++){
	236	Float_t r1=TMath::Sqrt(xc1[i]xc1[i]+yc1[i]yc1[i]);
	237	for(Int_t j=0;j<nrpL2;j++){
	238	Float_t diff = TMath::Abs(phi2[j]-phi1[i]);
	239	if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
	240	if(diff<fDiffPhiMax){
	241	Float_t r2=TMath::Sqrt(xc2[j]xc2[j]+yc2[j]yc2[j]);
	242	Float_t zr0=(r2zc1[i]-r1zc2[j])/(r2-r1);
	243	fZCombf->Fill(zr0);
	244	fZCombc->Fill(zr0);
	245	}
	246	}
	247	}
	248	delete [] xc1;
	249	delete [] yc1;
	250	delete [] zc1;
	251	delete [] phi1;
	252	delete [] xc2;
	253	delete [] yc2;
	254	delete [] zc2;
	255	delete [] phi2;
	256	if(fZCombc->GetEntries()==0){
	257	Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
	258	ResetHistograms();
	259	return fCurrentVertex;
	260	}
2257f27e	261	// else {
	262	// cout<<"Number of entries in hist. "<<fZCombc->GetEntries()<<endl;
	263	// }
0c6af5c9	264	Int_t bi = fZCombc->GetMaximumBin();
	265	Float_t centre = fZCombc->GetBinCenter(bi);
	266	Int_t n1 = static_cast<Int_t>((centre-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
	267	Int_t n2 = static_cast<Int_t>((centre+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
	268	Int_t niter = 0;
	269	Bool_t goon = kTRUE;
	270	Int_t num;
	271	while(goon){
	272	fZFound = 0.;
	273	fZsig = 0.;
	274	num=0;
	275	for(Int_t n=n1;n<=n2;n++){
	276	fZFound+=fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
	277	num+=static_cast<Int_t>(fZCombf->GetBinContent(n));
	278	fZsig+=fZCombf->GetBinCenter(n)fZCombf->GetBinCenter(n)fZCombf->GetBinContent(n);
	279	}
	280	if(num<2){
0c6af5c9	281	fZsig = 0.;
	282	}
	283	else {
979e3647	284	Float_t radi = fZsig/(num-1)-fZFound*fZFound/num/(num-1);
	285	if(radi>0.)fZsig=TMath::Sqrt(radi);
	286	else fZsig=0.;
0c6af5c9	287	fZFound/=num;
	288	}
	289	goon = TMath::Abs(TMath::Abs(fZFound-fZCombf->GetBinCenter(n1))-TMath::Abs(fZFound-fZCombf->GetBinCenter(n2)))>fTolerance;
0c6af5c9	290	n1 = static_cast<Int_t>((fZFound-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
	291	n2 = static_cast<Int_t>((fZFound+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
	292	niter++;
	293	if(niter>=10){
	294	goon = kFALSE;
	295	Warning("FindVertexForCurrentEvent","The procedure dows not converge\n");
	296	}
	297	}
2257f27e	298	// cout<<"Numer of Iterations "<<niter<<endl<<endl;
d681bb2d	299	fCurrentVertex = new AliESDVertex(fZFound,fZsig,num);
0c6af5c9	300	fCurrentVertex->SetTitle("vertexer: B");
	301	ResetHistograms();
	302	return fCurrentVertex;
	303	}
	304
	305	//_____________________________________________________________________
	306	void AliITSVertexerZ::ResetHistograms(){
	307	// delete TH1 data members
	308	if(fZCombc)delete fZCombc;
	309	if(fZCombf)delete fZCombf;
	310	fZCombc = 0;
	311	fZCombf = 0;
	312	}
	313
	314	//______________________________________________________________________
	315	void AliITSVertexerZ::FindVertices(){
	316	// computes the vertices of the events in the range FirstEvent - LastEvent
	317	AliRunLoader *rl = AliRunLoader::GetRunLoader();
	318	AliITSLoader* itsLoader = (AliITSLoader*) rl->GetLoader("ITSLoader");
	319	itsLoader->ReloadRecPoints();
	320	for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
	321	cout<<"Processing event "<<i<<endl;
	322	rl->GetEvent(i);
	323	FindVertexForCurrentEvent(i);
	324	if(fCurrentVertex){
	325	WriteCurrentVertex();
	326	}
	327	else {
	328	if(fDebug>0){
	329	cout<<"Vertex not found for event "<<i<<endl;
	330	cout<<"fZFound = "<<fZFound<<", fZsig= "<<fZsig<<endl;
	331	}
	332	}
	333	}
	334	}
	335
	336	//________________________________________________________
	337	void AliITSVertexerZ::PrintStatus() const {
	338	// Print current status
	339	cout <<"=======================================================\n";
	340	cout <<" First layer first and last modules: "<<fFirstL1<<", ";
	341	cout <<fLastL1<<endl;
	342	cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
	343	cout <<fLastL2<<endl;
	344	cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
	345	cout <<"Limits for Z histograms: "<<fLowLim<<"; "<<fHighLim<<endl;
	346	cout <<"Bin sizes for coarse and fine z histos "<<fStepCoarse<<"; "<<fStepFine<<endl;
	347	cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
	348	cout <<" Debug flag: "<<fDebug<<endl;
	349	cout <<"First event to be processed "<<fFirstEvent;
	350	cout <<"\n Last event to be processed "<<fLastEvent<<endl;
	351
	352
	353	cout <<"=======================================================\n";
	354	}
	355