[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<TBranch.h>
#include<TClonesArray.h>
#include<TH1.h>
#include <TString.h>
#include<TTree.h>
#include "AliESDVertex.h"
#include "AliITSgeomTGeo.h"
#include "AliITSRecPoint.h"
#include "AliITSRecPointContainer.h"
#include "AliITSZPoint.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
////////////////////////////////////////////////////////////////

using std::endl;
using std::cout;
ClassImp(AliITSVertexerZ)


//______________________________________________________________________
AliITSVertexerZ::AliITSVertexerZ():AliITSVertexer(),
fFirstL1(0),
fLastL1(0),
fFirstL2(0),
fLastL2(0),
fDiffPhiMax(0),
fZFound(0),
fZsig(0.),
fZCombc(0),
fLowLim(0.),
fHighLim(0.),
fStepCoarse(0),
fTolerance(0.),
fMaxIter(0),
fWindowWidth(0),
fSearchForPileup(kTRUE)
{
  // Default constructor
  SetDiffPhiMax();
  SetFirstLayerModules();
  SetSecondLayerModules();
  SetLowLimit();
  SetHighLimit();
  SetBinWidthCoarse();
  SetTolerance();
  SetPPsetting();
  ConfigIterations();
  SetWindowWidth();
}

//______________________________________________________________________
AliITSVertexerZ::AliITSVertexerZ(Float_t x0, Float_t y0):AliITSVertexer(),
fFirstL1(0),
fLastL1(0),
fFirstL2(0),
fLastL2(0),
fDiffPhiMax(0),
fZFound(0),
fZsig(0.),
fZCombc(0),
fLowLim(0.),
fHighLim(0.),
fStepCoarse(0),
fTolerance(0.),
fMaxIter(0),
fWindowWidth(0),
fSearchForPileup(kTRUE)
{
  // Standard Constructor
  SetDiffPhiMax();
  SetFirstLayerModules();
  SetSecondLayerModules();
  SetLowLimit();
  SetHighLimit();
  SetBinWidthCoarse();
  SetTolerance();
  SetPPsetting();
  ConfigIterations();
  SetWindowWidth();
  SetVtxStart((Double_t)x0,(Double_t)y0,0.);

}

//______________________________________________________________________
AliITSVertexerZ::~AliITSVertexerZ() {
  // Destructor
  delete fZCombc;
}

//______________________________________________________________________
void AliITSVertexerZ::ConfigIterations(Int_t noiter,Float_t *ptr){
  // configure the iterative procedure to gain efficiency for
  // pp events with very low multiplicity
  Float_t defaults[5]={0.02,0.05,0.1,0.2,0.3};
  fMaxIter=noiter;
  if(noiter>5){
    Error("ConfigIterations","Maximum number of iterations is 5\n");
    fMaxIter=5;
  }
  for(Int_t j=0;j<5;j++)fPhiDiffIter[j]=defaults[j];
  if(ptr)for(Int_t j=0;j<fMaxIter;j++)fPhiDiffIter[j]=ptr[j];
}

//______________________________________________________________________
Int_t AliITSVertexerZ::GetPeakRegion(TH1F*h, Int_t &binmin, Int_t &binmax){
  // Finds a region around a peak in the Z histogram
  // Case of 2 peaks is treated 
  Int_t imax=h->GetNbinsX();
  Float_t maxval=0;
  Int_t bi1=h->GetMaximumBin();
  Int_t bi2=0;
  for(Int_t i=imax;i>=1;i--){
    if(h->GetBinContent(i)>maxval){
      maxval=h->GetBinContent(i);
      bi2=i;
    }
  }
  Int_t npeaks=0;

  if(bi1==bi2){
    binmin=bi1-3;
    binmax=bi1+3;
    npeaks=1;
  }else{
    TH1F *copy = new TH1F(*h);
    copy->SetBinContent(bi1,0.);
    copy->SetBinContent(bi2,0.);
    Int_t l1=TMath::Max(bi1-3,1);
    Int_t l2=TMath::Min(bi1+3,h->GetNbinsX());
    Float_t cont1=copy->Integral(l1,l2);
    Int_t ll1=TMath::Max(bi2-3,1);
    Int_t ll2=TMath::Min(bi2+3,h->GetNbinsX());
    Float_t cont2=copy->Integral(ll1,ll2);
    if(cont1>cont2){
      binmin=l1;
      binmax=l2;
      npeaks=1;
    }
    if(cont2>cont1){
      binmin=ll1;
      binmax=ll2;
      npeaks=1;
    }
    if(cont1==cont2){
      binmin=l1;
      binmax=ll2;
      if(bi2-bi1==1) npeaks=1;
      else npeaks=2;
    }  
    delete copy;
  }    
  return npeaks;
}
//______________________________________________________________________
AliESDVertex* AliITSVertexerZ::FindVertexForCurrentEvent(TTree *itsClusterTree){
  // Defines the AliESDVertex for the current event
  VertexZFinder(itsClusterTree);
  Int_t ntrackl=0;
  for(Int_t iteraz=0;iteraz<fMaxIter;iteraz++){
    if(fCurrentVertex) ntrackl=fCurrentVertex->GetNContributors();
    if(!fCurrentVertex || ntrackl==0 || ntrackl==-1){
      Float_t diffPhiMaxOrig=fDiffPhiMax;
      fDiffPhiMax=GetPhiMaxIter(iteraz);
      VertexZFinder(itsClusterTree);
      fDiffPhiMax=diffPhiMaxOrig;
    }
  }
  if(fComputeMultiplicity) FindMultiplicity(itsClusterTree);
  return fCurrentVertex;
}  

//______________________________________________________________________
void AliITSVertexerZ::VertexZFinder(TTree *itsClusterTree){
  // Defines the AliESDVertex for the current event
  fCurrentVertex = 0;
  Double_t startPos[3]={GetNominalPos()[0],GetNominalPos()[1],GetNominalPos()[2]};
  Double_t startCov[6]={GetNominalCov()[0],GetNominalCov()[1],GetNominalCov()[2],
			GetNominalCov()[3],GetNominalCov()[4],GetNominalCov()[5]};
  ResetVertex();
  TClonesArray *itsRec  = 0;
  // lc1 and gc1 are local and global coordinates for layer 1
  Float_t gc1[3]={0.,0.,0.}; // ; for(Int_t ii=0; ii<3; ii++) gc1[ii]=0.;
  // lc2 and gc2 are local and global coordinates for layer 2
  Float_t gc2[3]={0.,0.,0.}; //; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
  AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
  rpcont->FetchClusters(0,itsClusterTree);
  if(!rpcont->IsSPDActive()){
    AliWarning("Null pointer for RecPoints branch, vertex not calculated");
    ResetHistograms();
    fCurrentVertex = new AliESDVertex(startPos,startCov,99999.,-2);
    return;
  }

  Int_t nrpL1 = 0;
  Int_t nrpL2 = 0;
  nrpL1=rpcont->GetNClustersInLayerFast(1);
  nrpL2=rpcont->GetNClustersInLayerFast(2);

  if(nrpL1 == 0 || nrpL2 == 0){
    AliDebug(1,Form("No RecPoints in at least one SPD layer (%d %d)",nrpL1,nrpL2));
    ResetHistograms();
    fCurrentVertex = new AliESDVertex(startPos,startCov,99999.,-2);
    return;
  }
  // Force a coarse bin size of 200 microns if the number of clusters on layer 2
  // is low
  if(nrpL2<fPPsetting[0])SetBinWidthCoarse(fPPsetting[1]);
  // By default nbincoarse=(10+10)/0.01=2000
  Int_t nbincoarse = static_cast<Int_t>((fHighLim-fLowLim)/fStepCoarse);
  if(fZCombc)delete fZCombc;
  fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);

 /* Test the ffect of mutiple scatternig on error. Negligible
  // Multiple scattering
  Float_t beta=1.,pmed=0.875; //pmed=875 MeV (for tracks with dphi<0.01 rad)
  Float_t beta2=beta*beta;
  Float_t p2=pmed*pmed;
  Float_t rBP=3; //Beam Pipe radius = 3cm
  Float_t dBP=0.08/35.3; // 800 um of Be
  Float_t dL1=0.01; //approx. 1% of radiation length  
  Float_t theta2BP=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(dBP);
  Float_t theta2L1=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(dL1);
*/
  Int_t nEntriesMod[kNSPDMod];
  TClonesArray* recpArr[kNSPDMod];
  for(Int_t modul=0; modul<kNSPDMod; ++modul) {
    if(!fUseModule[modul]) {
      nEntriesMod[modul]=0;
      recpArr[modul]=0;
    } else {
      recpArr[modul]=rpcont->UncheckedGetClusters(modul);
      nEntriesMod[modul]=recpArr[modul]->GetEntriesFast();
    }
  }
      
  Int_t maxdim=TMath::Min(nrpL1*nrpL2,50000);  // temporary; to limit the size in PbPb
  static TClonesArray points("AliITSZPoint",maxdim);
  Int_t nopoints =0;
  for(Int_t modul1= fFirstL1; modul1<=fLastL1;modul1++){   // Loop on modules of layer 1
    if(!fUseModule[modul1]) continue;
    UShort_t ladder=int(modul1/4)+1;  // ladders are numbered starting from 1
    TClonesArray *prpl1=recpArr[modul1]; //rpcont->UncheckedGetClusters(modul1);
    Int_t nrecp1 = nEntriesMod[modul1]; //prpl1->GetEntries();
    for(Int_t j1=0;j1<nrecp1;j1++){
      AliITSRecPoint *recp1 = (AliITSRecPoint*)prpl1->At(j1);
      recp1->GetGlobalXYZ(gc1);
      gc1[0]-=GetNominalPos()[0]; // Possible beam offset in the bending plane
      gc1[1]-=GetNominalPos()[1]; //   "               "
      Float_t phi1 = TMath::ATan2(gc1[1],gc1[0]);
      if(phi1<0)phi1+=TMath::TwoPi();
      for(Int_t ladl2=0 ; ladl2<fLadOnLay2*2+1;ladl2++){
	for(Int_t k=0;k<4;k++){
 	  Int_t ladmod=fLadders[ladder-1]+ladl2;
	  if(ladmod>AliITSgeomTGeo::GetNLadders(2)) ladmod=ladmod-AliITSgeomTGeo::GetNLadders(2);
	  Int_t modul2=AliITSgeomTGeo::GetModuleIndex(2,ladmod,k+1);
	  if(modul2<0)continue;
	  if(!fUseModule[modul2]) continue;
	  itsRec=recpArr[modul2]; // rpcont->UncheckedGetClusters(modul2);
	  Int_t nrecp2 = nEntriesMod[modul2]; // itsRec->GetEntries();
	  for(Int_t j2=0;j2<nrecp2;j2++){
	    AliITSRecPoint *recp2 = (AliITSRecPoint*)itsRec->At(j2);
	    recp2->GetGlobalXYZ(gc2);
	    gc2[0]-=GetNominalPos()[0];
	    gc2[1]-=GetNominalPos()[1];
	    Float_t phi2 = TMath::ATan2(gc2[1],gc2[0]);
	    if(phi2<0)phi2+=TMath::TwoPi();

	    Float_t diff = TMath::Abs(phi2-phi1); 
	    if(diff>TMath::Pi())diff=TMath::TwoPi()-diff;
	    if(diff<fDiffPhiMax){
	      Float_t r1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]);
	      Float_t zc1=gc1[2];
	      Float_t erz1=recp1->GetSigmaZ2();
	      Float_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]);
	      Float_t zc2=gc2[2];
	      Float_t erz2=recp2->GetSigmaZ2();
	      //	Float_t tgth=(zc2[j]-zc1[i])/(r2-r1); // slope (used for multiple scattering)
	      Float_t zr0=(r2*zc1-r1*zc2)/(r2-r1); //Z @ null radius
	      Float_t ezr0q=(r2*r2*erz1+r1*r1*erz2)/((r2-r1)*(r2-r1)); //error on Z @ null radius
	      /*
	      // Multiple scattering
	      ezr0q+=r1*r1*(1+tgth*tgth)*theta2L1/2; // multiple scattering in layer 1
	      ezr0q+=rBP*rBP*(1+tgth*tgth)*theta2BP/2; // multiple scattering in beam pipe
	      */
	      if(nopoints<maxdim) new(points[nopoints++])AliITSZPoint(zr0,ezr0q);	      
	      fZCombc->Fill(zr0);
	    }
	  }
	}
      }
    }
  }

  points.Sort();

  Double_t contents = fZCombc->GetEntries()- fZCombc->GetBinContent(0)-fZCombc->GetBinContent(nbincoarse+1);
  if(contents<1.){
    //    Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
    ResetHistograms();
    fCurrentVertex = new AliESDVertex(startPos,startCov,99999.,-1);
    points.Clear();
    return;
  }

  TH1F *hc = fZCombc;

  
  if(hc->GetBinContent(hc->GetMaximumBin())<3)hc->Rebin(4);
  Int_t binmin,binmax;
  Int_t nPeaks=GetPeakRegion(hc,binmin,binmax);   
  if(nPeaks==2)AliDebug(2,"2 peaks found");
  Float_t zm =0.;
  Float_t ezm =0.;
  Float_t lim1 = hc->GetBinLowEdge(binmin);
  Float_t lim2 = hc->GetBinLowEdge(binmax)+hc->GetBinWidth(binmax);
  Float_t widthSR=lim2-lim1;

  if(nPeaks ==1 && (lim2-lim1)<fWindowWidth){
    Float_t c=(lim1+lim2)/2.;
    lim1=c-fWindowWidth/2.;
    lim2=c+fWindowWidth/2.;
  }
  Int_t niter = 0, ncontr=0;
  do {
    // symmetrization
    if(zm  !=0.){
      Float_t semilarg=TMath::Min((lim2-zm),(zm-lim1));
      lim1=zm - semilarg;
      lim2=zm + semilarg;
    }

    zm=0.;
    ezm=0.;
    ncontr=0;
    for(Int_t i =0; i<points.GetEntriesFast(); i++){
      AliITSZPoint* p=(AliITSZPoint*)points.UncheckedAt(i);
      if(p->GetZ()>lim1 && p->GetZ()<lim2){
        Float_t deno = p->GetErrZ();
        zm+=p->GetZ()/deno;
        ezm+=1./deno;
        ncontr++;
      }
    }
    if(ezm>0) {
      zm/=ezm;
      ezm=TMath::Sqrt(1./ezm);
    }
    niter++;
  } while(niter<10 && TMath::Abs((zm-lim1)-(lim2-zm))>fTolerance);
  if(nPeaks==2) ezm=widthSR;
  Double_t position[3]={GetNominalPos()[0],GetNominalPos()[1],zm};
  Double_t covmatrix[6]={GetNominalCov()[0],0.,GetNominalCov()[2],0.,0.,ezm};
  fCurrentVertex = new AliESDVertex(position,covmatrix,99999.,ncontr);
  fCurrentVertex->SetTitle("vertexer: Z");
  fCurrentVertex->SetDispersion(fDiffPhiMax);
  fNoVertices=1;
  points.Clear();
  if(fSearchForPileup && ncontr>fMinTrackletsForPilup){ 
    Float_t secPeakPos;
    Int_t ncontr2=FindSecondPeak(fZCombc,binmin,binmax,secPeakPos);
    if(ncontr2>=fMinTrackletsForPilup){ 
      fIsPileup=kTRUE;
      fNoVertices=2;
      fZpuv=secPeakPos;
      fNTrpuv=ncontr2;
      AliESDVertex secondVert(secPeakPos,0.1,ncontr2);
      fVertArray = new AliESDVertex[2];
      fVertArray[0]=(*fCurrentVertex);
      fVertArray[1]=secondVert;
    }
  }
  if(fNoVertices==1){
    fVertArray = new AliESDVertex[1];
    fVertArray[0]=(*fCurrentVertex);	  
  }
  
  ResetHistograms();
  return;
}

//_____________________________________________________________________
Int_t AliITSVertexerZ::FindSecondPeak(TH1F* h, Int_t binmin,Int_t binmax, Float_t& secPeakPos){ 
  // Resets bin contents between binmin and binmax and then search 
  // for a second peak position 
  for(Int_t i=binmin-1;i<=binmax+1;i++){
    h->SetBinContent(i,0.);
  }
  Int_t secPeakBin=h->GetMaximumBin();
  secPeakPos=h->GetBinCenter(secPeakBin);
  Int_t secPeakCont=(Int_t)h->GetBinContent(secPeakBin);
  secPeakCont+=(Int_t)h->GetBinContent(secPeakBin-1);
  secPeakCont+=(Int_t)h->GetBinContent(secPeakBin+1);  
  secPeakCont+=(Int_t)h->GetBinContent(secPeakBin-2);
  secPeakCont+=(Int_t)h->GetBinContent(secPeakBin+2);  
  return secPeakCont;
}

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

//________________________________________________________
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 z histos "<<fStepCoarse<<endl;
  cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
  if(fZCombc){
    cout<<"fZCombc exists - entries="<<fZCombc->GetEntries()<<endl;
  }
  else{
    cout<<"fZCombc does not exist\n";
  }
 
  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	**************************************************************************/
d531b7a1	15	#include "AliITSVertexerZ.h"
0c6af5c9	16	#include<TBranch.h>
7299b7a8	17	#include<TClonesArray.h>
0c6af5c9	18	#include<TH1.h>
d531b7a1	19	#include <TString.h>
0c6af5c9	20	#include<TTree.h>
308c2f7c	21	#include "AliESDVertex.h"
32e63e47	22	#include "AliITSgeomTGeo.h"
d531b7a1	23	#include "AliITSRecPoint.h"
b21c1af0	24	#include "AliITSRecPointContainer.h"
7299b7a8	25	#include "AliITSZPoint.h"
0c6af5c9	26
	27	/////////////////////////////////////////////////////////////////
	28	// this class implements a fast method to determine
	29	// the Z coordinate of the primary vertex
	30	// for p-p collisions it seems to give comparable or better results
	31	// with respect to what obtained with AliITSVertexerPPZ
	32	// It can be used successfully with Pb-Pb collisions
	33	////////////////////////////////////////////////////////////////
	34
fe7d86eb	35	using std::endl;
fe7d86eb	36	using std::cout;
0c6af5c9	37	ClassImp(AliITSVertexerZ)
	38
	39
	40
	41	//______________________________________________________________________
8221b41b	42	AliITSVertexerZ::AliITSVertexerZ():AliITSVertexer(),
	43	fFirstL1(0),
	44	fLastL1(0),
	45	fFirstL2(0),
	46	fLastL2(0),
	47	fDiffPhiMax(0),
8221b41b	48	fZFound(0),
	49	fZsig(0.),
	50	fZCombc(0),
8221b41b	51	fLowLim(0.),
	52	fHighLim(0.),
	53	fStepCoarse(0),
8221b41b	54	fTolerance(0.),
7299b7a8	55	fMaxIter(0),
1e6316d8	56	fWindowWidth(0),
	57	fSearchForPileup(kTRUE)
	58	{
0c6af5c9	59	// Default constructor
7299b7a8	60	SetDiffPhiMax();
	61	SetFirstLayerModules();
	62	SetSecondLayerModules();
	63	SetLowLimit();
	64	SetHighLimit();
	65	SetBinWidthCoarse();
	66	SetTolerance();
	67	SetPPsetting();
6fd990e3	68	ConfigIterations();
7299b7a8	69	SetWindowWidth();
0c6af5c9	70	}
	71
	72	//______________________________________________________________________
308c2f7c	73	AliITSVertexerZ::AliITSVertexerZ(Float_t x0, Float_t y0):AliITSVertexer(),
8221b41b	74	fFirstL1(0),
	75	fLastL1(0),
	76	fFirstL2(0),
	77	fLastL2(0),
	78	fDiffPhiMax(0),
8221b41b	79	fZFound(0),
	80	fZsig(0.),
	81	fZCombc(0),
8221b41b	82	fLowLim(0.),
	83	fHighLim(0.),
	84	fStepCoarse(0),
8221b41b	85	fTolerance(0.),
7299b7a8	86	fMaxIter(0),
1e6316d8	87	fWindowWidth(0),
	88	fSearchForPileup(kTRUE)
	89	{
0c6af5c9	90	// Standard Constructor
0c6af5c9	91	SetDiffPhiMax();
0c6af5c9	92	SetFirstLayerModules();
0c6af5c9	93	SetSecondLayerModules();
0c6af5c9	94	SetLowLimit();
	95	SetHighLimit();
	96	SetBinWidthCoarse();
0c6af5c9	97	SetTolerance();
ecc64c3f	98	SetPPsetting();
6fd990e3	99	ConfigIterations();
7299b7a8	100	SetWindowWidth();
60b9526b	101	SetVtxStart((Double_t)x0,(Double_t)y0,0.);
0c6af5c9	102
	103	}
	104
0c6af5c9	105	//______________________________________________________________________
0c6af5c9	106	AliITSVertexerZ::~AliITSVertexerZ() {
ecc64c3f	107	// Destructor
ecc64c3f	108	delete fZCombc;
0c6af5c9	109	}
0c6af5c9	110
6fd990e3	111	//______________________________________________________________________
	112	void AliITSVertexerZ::ConfigIterations(Int_t noiter,Float_t *ptr){
	113	// configure the iterative procedure to gain efficiency for
	114	// pp events with very low multiplicity
bf2e0ad4	115	Float_t defaults[5]={0.02,0.05,0.1,0.2,0.3};
6fd990e3	116	fMaxIter=noiter;
	117	if(noiter>5){
	118	Error("ConfigIterations","Maximum number of iterations is 5\n");
	119	fMaxIter=5;
	120	}
	121	for(Int_t j=0;j<5;j++)fPhiDiffIter[j]=defaults[j];
	122	if(ptr)for(Int_t j=0;j<fMaxIter;j++)fPhiDiffIter[j]=ptr[j];
	123	}
	124
debb3dde	125	//______________________________________________________________________
8c32ba44	126	Int_t AliITSVertexerZ::GetPeakRegion(TH1F*h, Int_t &binmin, Int_t &binmax){
debb3dde	127	// Finds a region around a peak in the Z histogram
	128	// Case of 2 peaks is treated
	129	Int_t imax=h->GetNbinsX();
	130	Float_t maxval=0;
	131	Int_t bi1=h->GetMaximumBin();
	132	Int_t bi2=0;
	133	for(Int_t i=imax;i>=1;i--){
	134	if(h->GetBinContent(i)>maxval){
	135	maxval=h->GetBinContent(i);
	136	bi2=i;
	137	}
	138	}
	139	Int_t npeaks=0;
	140
	141	if(bi1==bi2){
	142	binmin=bi1-3;
	143	binmax=bi1+3;
	144	npeaks=1;
	145	}else{
	146	TH1F copy = new TH1F(h);
	147	copy->SetBinContent(bi1,0.);
	148	copy->SetBinContent(bi2,0.);
	149	Int_t l1=TMath::Max(bi1-3,1);
	150	Int_t l2=TMath::Min(bi1+3,h->GetNbinsX());
	151	Float_t cont1=copy->Integral(l1,l2);
	152	Int_t ll1=TMath::Max(bi2-3,1);
	153	Int_t ll2=TMath::Min(bi2+3,h->GetNbinsX());
	154	Float_t cont2=copy->Integral(ll1,ll2);
	155	if(cont1>cont2){
	156	binmin=l1;
	157	binmax=l2;
	158	npeaks=1;
	159	}
	160	if(cont2>cont1){
	161	binmin=ll1;
	162	binmax=ll2;
	163	npeaks=1;
	164	}
	165	if(cont1==cont2){
	166	binmin=l1;
	167	binmax=ll2;
	168	if(bi2-bi1==1) npeaks=1;
	169	else npeaks=2;
	170	}
	171	delete copy;
	172	}
	173	return npeaks;
	174	}
0c6af5c9	175	//______________________________________________________________________
308c2f7c	176	AliESDVertex* AliITSVertexerZ::FindVertexForCurrentEvent(TTree *itsClusterTree){
d681bb2d	177	// Defines the AliESDVertex for the current event
308c2f7c	178	VertexZFinder(itsClusterTree);
6fd990e3	179	Int_t ntrackl=0;
	180	for(Int_t iteraz=0;iteraz<fMaxIter;iteraz++){
	181	if(fCurrentVertex) ntrackl=fCurrentVertex->GetNContributors();
	182	if(!fCurrentVertex \|\| ntrackl==0 \|\| ntrackl==-1){
	183	Float_t diffPhiMaxOrig=fDiffPhiMax;
	184	fDiffPhiMax=GetPhiMaxIter(iteraz);
308c2f7c	185	VertexZFinder(itsClusterTree);
6fd990e3	186	fDiffPhiMax=diffPhiMaxOrig;
	187	}
	188	}
1ff24d0a	189	if(fComputeMultiplicity) FindMultiplicity(itsClusterTree);
6fd990e3	190	return fCurrentVertex;
	191	}
	192
	193	//______________________________________________________________________
308c2f7c	194	void AliITSVertexerZ::VertexZFinder(TTree *itsClusterTree){
6fd990e3	195	// Defines the AliESDVertex for the current event
0c6af5c9	196	fCurrentVertex = 0;
5279841d	197	Double_t startPos[3]={GetNominalPos()[0],GetNominalPos()[1],GetNominalPos()[2]};
	198	Double_t startCov[6]={GetNominalCov()[0],GetNominalCov()[1],GetNominalCov()[2],
	199	GetNominalCov()[3],GetNominalCov()[4],GetNominalCov()[5]};
6b4d9537	200	ResetVertex();
0c6af5c9	201	TClonesArray *itsRec = 0;
27167524	202	// lc1 and gc1 are local and global coordinates for layer 1
54f95b99	203	Float_t gc1[3]={0.,0.,0.}; // ; for(Int_t ii=0; ii<3; ii++) gc1[ii]=0.;
ecc64c3f	204	// lc2 and gc2 are local and global coordinates for layer 2
54f95b99	205	Float_t gc2[3]={0.,0.,0.}; //; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
b21c1af0	206	AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
3457dd2d	207	rpcont->FetchClusters(0,itsClusterTree);
b21c1af0	208	if(!rpcont->IsSPDActive()){
fea0ad94	209	AliWarning("Null pointer for RecPoints branch, vertex not calculated");
fea0ad94	210	ResetHistograms();
5279841d	211	fCurrentVertex = new AliESDVertex(startPos,startCov,99999.,-2);
fea0ad94	212	return;
fea0ad94	213	}
0c6af5c9	214
0c6af5c9	215	Int_t nrpL1 = 0;
0c6af5c9	216	Int_t nrpL2 = 0;
1cc75a0b	217	nrpL1=rpcont->GetNClustersInLayerFast(1);
1cc75a0b	218	nrpL2=rpcont->GetNClustersInLayerFast(2);
27167524	219
0c6af5c9	220	if(nrpL1 == 0 \|\| nrpL2 == 0){
fea0ad94	221	AliDebug(1,Form("No RecPoints in at least one SPD layer (%d %d)",nrpL1,nrpL2));
0c6af5c9	222	ResetHistograms();
5279841d	223	fCurrentVertex = new AliESDVertex(startPos,startCov,99999.,-2);
6fd990e3	224	return;
0c6af5c9	225	}
ecc64c3f	226	// Force a coarse bin size of 200 microns if the number of clusters on layer 2
	227	// is low
	228	if(nrpL2<fPPsetting[0])SetBinWidthCoarse(fPPsetting[1]);
ecc64c3f	229	// By default nbincoarse=(10+10)/0.01=2000
ecc64c3f	230	Int_t nbincoarse = static_cast<Int_t>((fHighLim-fLowLim)/fStepCoarse);
ecc64c3f	231	if(fZCombc)delete fZCombc;
ecc64c3f	232	fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);
7299b7a8	233
27167524	234	/* Test the ffect of mutiple scatternig on error. Negligible
7299b7a8	235	// Multiple scattering
	236	Float_t beta=1.,pmed=0.875; //pmed=875 MeV (for tracks with dphi<0.01 rad)
	237	Float_t beta2=beta*beta;
	238	Float_t p2=pmed*pmed;
	239	Float_t rBP=3; //Beam Pipe radius = 3cm
	240	Float_t dBP=0.08/35.3; // 800 um of Be
	241	Float_t dL1=0.01; //approx. 1% of radiation length
	242	Float_t theta2BP=14.114.1/(beta2p21e6)TMath::Abs(dBP);
	243	Float_t theta2L1=14.114.1/(beta2p21e6)TMath::Abs(dL1);
	244	*/
54f95b99	245	Int_t nEntriesMod[kNSPDMod];
	246	TClonesArray* recpArr[kNSPDMod];
	247	for(Int_t modul=0; modul<kNSPDMod; ++modul) {
	248	if(!fUseModule[modul]) {
	249	nEntriesMod[modul]=0;
	250	recpArr[modul]=0;
	251	} else {
	252	recpArr[modul]=rpcont->UncheckedGetClusters(modul);
	253	nEntriesMod[modul]=recpArr[modul]->GetEntriesFast();
	254	}
	255	}
	256
27167524	257	Int_t maxdim=TMath::Min(nrpL1*nrpL2,50000); // temporary; to limit the size in PbPb
f5f6da22	258	static TClonesArray points("AliITSZPoint",maxdim);
7299b7a8	259	Int_t nopoints =0;
27167524	260	for(Int_t modul1= fFirstL1; modul1<=fLastL1;modul1++){ // Loop on modules of layer 1
8c42830a	261	if(!fUseModule[modul1]) continue;
27167524	262	UShort_t ladder=int(modul1/4)+1; // ladders are numbered starting from 1
54f95b99	263	TClonesArray *prpl1=recpArr[modul1]; //rpcont->UncheckedGetClusters(modul1);
54f95b99	264	Int_t nrecp1 = nEntriesMod[modul1]; //prpl1->GetEntries();
27167524	265	for(Int_t j1=0;j1<nrecp1;j1++){
54f95b99	266	AliITSRecPoint recp1 = (AliITSRecPoint)prpl1->At(j1);
54f95b99	267	recp1->GetGlobalXYZ(gc1);
308c2f7c	268	gc1[0]-=GetNominalPos()[0]; // Possible beam offset in the bending plane
308c2f7c	269	gc1[1]-=GetNominalPos()[1]; // " "
27167524	270	Float_t phi1 = TMath::ATan2(gc1[1],gc1[0]);
54f95b99	271	if(phi1<0)phi1+=TMath::TwoPi();
27167524	272	for(Int_t ladl2=0 ; ladl2<fLadOnLay2*2+1;ladl2++){
	273	for(Int_t k=0;k<4;k++){
	274	Int_t ladmod=fLadders[ladder-1]+ladl2;
32e63e47	275	if(ladmod>AliITSgeomTGeo::GetNLadders(2)) ladmod=ladmod-AliITSgeomTGeo::GetNLadders(2);
32e63e47	276	Int_t modul2=AliITSgeomTGeo::GetModuleIndex(2,ladmod,k+1);
e917bef4	277	if(modul2<0)continue;
8c42830a	278	if(!fUseModule[modul2]) continue;
54f95b99	279	itsRec=recpArr[modul2]; // rpcont->UncheckedGetClusters(modul2);
54f95b99	280	Int_t nrecp2 = nEntriesMod[modul2]; // itsRec->GetEntries();
27167524	281	for(Int_t j2=0;j2<nrecp2;j2++){
54f95b99	282	AliITSRecPoint recp2 = (AliITSRecPoint)itsRec->At(j2);
54f95b99	283	recp2->GetGlobalXYZ(gc2);
308c2f7c	284	gc2[0]-=GetNominalPos()[0];
308c2f7c	285	gc2[1]-=GetNominalPos()[1];
27167524	286	Float_t phi2 = TMath::ATan2(gc2[1],gc2[0]);
54f95b99	287	if(phi2<0)phi2+=TMath::TwoPi();
27167524	288
27167524	289	Float_t diff = TMath::Abs(phi2-phi1);
54f95b99	290	if(diff>TMath::Pi())diff=TMath::TwoPi()-diff;
27167524	291	if(diff<fDiffPhiMax){
54f95b99	292	Float_t r1=TMath::Sqrt(gc1[0]gc1[0]+gc1[1]gc1[1]);
	293	Float_t zc1=gc1[2];
	294	Float_t erz1=recp1->GetSigmaZ2();
	295	Float_t r2=TMath::Sqrt(gc2[0]gc2[0]+gc2[1]gc2[1]);
	296	Float_t zc2=gc2[2];
	297	Float_t erz2=recp2->GetSigmaZ2();
27167524	298	// Float_t tgth=(zc2[j]-zc1[i])/(r2-r1); // slope (used for multiple scattering)
27167524	299	Float_t zr0=(r2zc1-r1zc2)/(r2-r1); //Z @ null radius
54f95b99	300	Float_t ezr0q=(r2r2erz1+r1r1erz2)/((r2-r1)*(r2-r1)); //error on Z @ null radius
	301	/*
	302	// Multiple scattering
	303	ezr0q+=r1r1(1+tgthtgth)theta2L1/2; // multiple scattering in layer 1
	304	ezr0q+=rBPrBP(1+tgthtgth)theta2BP/2; // multiple scattering in beam pipe
	305	*/
f5f6da22	306	if(nopoints<maxdim) new(points[nopoints++])AliITSZPoint(zr0,ezr0q);
27167524	307	fZCombc->Fill(zr0);
	308	}
	309	}
27167524	310	}
0c6af5c9	311	}
	312	}
	313	}
7299b7a8	314
f5f6da22	315	points.Sort();
7299b7a8	316
6fd990e3	317	Double_t contents = fZCombc->GetEntries()- fZCombc->GetBinContent(0)-fZCombc->GetBinContent(nbincoarse+1);
	318	if(contents<1.){
	319	// Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
0c6af5c9	320	ResetHistograms();
5279841d	321	fCurrentVertex = new AliESDVertex(startPos,startCov,99999.,-1);
f5f6da22	322	points.Clear();
6fd990e3	323	return;
0c6af5c9	324	}
ecc64c3f	325
ecc64c3f	326	TH1F *hc = fZCombc;
ecc64c3f	327
ecc64c3f	328
8c42830a	329	if(hc->GetBinContent(hc->GetMaximumBin())<3)hc->Rebin(4);
debb3dde	330	Int_t binmin,binmax;
debb3dde	331	Int_t nPeaks=GetPeakRegion(hc,binmin,binmax);
4be14da5	332	if(nPeaks==2)AliDebug(2,"2 peaks found");
7299b7a8	333	Float_t zm =0.;
	334	Float_t ezm =0.;
	335	Float_t lim1 = hc->GetBinLowEdge(binmin);
	336	Float_t lim2 = hc->GetBinLowEdge(binmax)+hc->GetBinWidth(binmax);
4be14da5	337	Float_t widthSR=lim2-lim1;
7299b7a8	338
	339	if(nPeaks ==1 && (lim2-lim1)<fWindowWidth){
	340	Float_t c=(lim1+lim2)/2.;
	341	lim1=c-fWindowWidth/2.;
	342	lim2=c+fWindowWidth/2.;
ecc64c3f	343	}
7299b7a8	344	Int_t niter = 0, ncontr=0;
	345	do {
	346	// symmetrization
	347	if(zm !=0.){
	348	Float_t semilarg=TMath::Min((lim2-zm),(zm-lim1));
	349	lim1=zm - semilarg;
	350	lim2=zm + semilarg;
0c6af5c9	351	}
7299b7a8	352
	353	zm=0.;
	354	ezm=0.;
	355	ncontr=0;
896141ab	356	for(Int_t i =0; i<points.GetEntriesFast(); i++){
f5f6da22	357	AliITSZPoint* p=(AliITSZPoint*)points.UncheckedAt(i);
7299b7a8	358	if(p->GetZ()>lim1 && p->GetZ()<lim2){
	359	Float_t deno = p->GetErrZ();
	360	zm+=p->GetZ()/deno;
	361	ezm+=1./deno;
	362	ncontr++;
ecc64c3f	363	}
ecc64c3f	364	}
2b8b4286	365	if(ezm>0) {
	366	zm/=ezm;
	367	ezm=TMath::Sqrt(1./ezm);
	368	}
7299b7a8	369	niter++;
7299b7a8	370	} while(niter<10 && TMath::Abs((zm-lim1)-(lim2-zm))>fTolerance);
4be14da5	371	if(nPeaks==2) ezm=widthSR;
5279841d	372	Double_t position[3]={GetNominalPos()[0],GetNominalPos()[1],zm};
	373	Double_t covmatrix[6]={GetNominalCov()[0],0.,GetNominalCov()[2],0.,0.,ezm};
	374	fCurrentVertex = new AliESDVertex(position,covmatrix,99999.,ncontr);
12e3ead8	375	fCurrentVertex->SetTitle("vertexer: Z");
bf2e0ad4	376	fCurrentVertex->SetDispersion(fDiffPhiMax);
6b4d9537	377	fNoVertices=1;
f5f6da22	378	points.Clear();
1e6316d8	379	if(fSearchForPileup && ncontr>fMinTrackletsForPilup){
8c32ba44	380	Float_t secPeakPos;
	381	Int_t ncontr2=FindSecondPeak(fZCombc,binmin,binmax,secPeakPos);
	382	if(ncontr2>=fMinTrackletsForPilup){
	383	fIsPileup=kTRUE;
6b4d9537	384	fNoVertices=2;
8c32ba44	385	fZpuv=secPeakPos;
8c32ba44	386	fNTrpuv=ncontr2;
6b4d9537	387	AliESDVertex secondVert(secPeakPos,0.1,ncontr2);
	388	fVertArray = new AliESDVertex[2];
	389	fVertArray[0]=(*fCurrentVertex);
	390	fVertArray[1]=secondVert;
8c32ba44	391	}
2c3a7cd0	392	}
6b4d9537	393	if(fNoVertices==1){
	394	fVertArray = new AliESDVertex[1];
	395	fVertArray[0]=(*fCurrentVertex);
	396	}
	397
0c6af5c9	398	ResetHistograms();
6fd990e3	399	return;
0c6af5c9	400	}
0c6af5c9	401
8c32ba44	402	//_____________________________________________________________________
f510fd70	403	Int_t AliITSVertexerZ::FindSecondPeak(TH1F* h, Int_t binmin,Int_t binmax, Float_t& secPeakPos){
	404	// Resets bin contents between binmin and binmax and then search
	405	// for a second peak position
8c32ba44	406	for(Int_t i=binmin-1;i<=binmax+1;i++){
	407	h->SetBinContent(i,0.);
	408	}
	409	Int_t secPeakBin=h->GetMaximumBin();
	410	secPeakPos=h->GetBinCenter(secPeakBin);
70a3283c	411	Int_t secPeakCont=(Int_t)h->GetBinContent(secPeakBin);
	412	secPeakCont+=(Int_t)h->GetBinContent(secPeakBin-1);
	413	secPeakCont+=(Int_t)h->GetBinContent(secPeakBin+1);
	414	secPeakCont+=(Int_t)h->GetBinContent(secPeakBin-2);
	415	secPeakCont+=(Int_t)h->GetBinContent(secPeakBin+2);
8c32ba44	416	return secPeakCont;
	417	}
	418
0c6af5c9	419	//_____________________________________________________________________
	420	void AliITSVertexerZ::ResetHistograms(){
	421	// delete TH1 data members
	422	if(fZCombc)delete fZCombc;
0c6af5c9	423	fZCombc = 0;
0c6af5c9	424	}
0c6af5c9	425
0c6af5c9	426	//________________________________________________________
	427	void AliITSVertexerZ::PrintStatus() const {
	428	// Print current status
	429	cout <<"=======================================================\n";
	430	cout <<" First layer first and last modules: "<<fFirstL1<<", ";
	431	cout <<fLastL1<<endl;
	432	cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
	433	cout <<fLastL2<<endl;
	434	cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
	435	cout <<"Limits for Z histograms: "<<fLowLim<<"; "<<fHighLim<<endl;
7299b7a8	436	cout <<"Bin sizes for coarse z histos "<<fStepCoarse<<endl;
0c6af5c9	437	cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
ecc64c3f	438	if(fZCombc){
	439	cout<<"fZCombc exists - entries="<<fZCombc->GetEntries()<<endl;
	440	}
	441	else{
	442	cout<<"fZCombc does not exist\n";
	443	}
0c6af5c9	444
	445	cout <<"=======================================================\n";
	446	}
	447