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