[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
////////////////////////////////////////////////////////////////

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) {
  // 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) {
  // 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;
    }
  }
  FindMultiplicity(itsClusterTree);
  return fCurrentVertex;
}  

//______________________________________________________________________
void AliITSVertexerZ::VertexZFinder(TTree *itsClusterTree){
  // Defines the AliESDVertex for the current event
  fCurrentVertex = 0;
  ResetVertex();
  TClonesArray *itsRec  = 0;
  // lc1 and gc1 are local and global coordinates for layer 1
  Float_t lc1[3]; for(Int_t ii=0; ii<3; ii++) lc1[ii]=0.;
  Float_t gc1[3]; for(Int_t ii=0; ii<3; ii++) gc1[ii]=0.;
  // lc2 and gc2 are local and global coordinates for layer 2
  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.;
  AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
  itsRec=rpcont->FetchClusters(0,itsClusterTree);
  if(!rpcont->IsSPDActive()){
    AliWarning("Null pointer for RecPoints branch, vertex not calculated");
    ResetHistograms();
    fCurrentVertex = new AliESDVertex(0.,5.3,-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(0.,5.3,-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 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=rpcont->UncheckedGetClusters(modul1);
    Int_t nrecp1 = prpl1->GetEntries();
    for(Int_t j1=0;j1<nrecp1;j1++){
      AliITSRecPoint *recp = (AliITSRecPoint*)prpl1->At(j1);
      recp->GetGlobalXYZ(gc1);
      gc1[0]-=GetNominalPos()[0]; // Possible beam offset in the bending plane
      gc1[1]-=GetNominalPos()[1]; //   "               "
      Float_t r1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]);
      Float_t phi1 = TMath::ATan2(gc1[1],gc1[0]);
      if(phi1<0)phi1+=2*TMath::Pi();
      Float_t zc1=gc1[2];
      Float_t erz1=recp->GetSigmaZ2();
      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(!fUseModule[modul2]) continue;
	  itsRec=rpcont->UncheckedGetClusters(modul2);
	  Int_t nrecp2 = itsRec->GetEntries();
	  for(Int_t j2=0;j2<nrecp2;j2++){
	    recp = (AliITSRecPoint*)itsRec->At(j2);
	    recp->GetGlobalXYZ(gc2);
	    gc2[0]-=GetNominalPos()[0];
	    gc2[1]-=GetNominalPos()[1];
	    Float_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]);
	    Float_t phi2 = TMath::ATan2(gc2[1],gc2[0]);
	    if(phi2<0)phi2+=2*TMath::Pi();
	    Float_t zc2=gc2[2];
	    Float_t erz2=recp->GetSigmaZ2();

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