[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;
    }
  }
  if(fComputeMultiplicity) 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]={0.,0.,0.}; // for(Int_t ii=0; ii<3; ii++) lc1[ii]=0.;
  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 lc2[3]={0.,0.,0.}; // ; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
  Float_t gc2[3]={0.,0.,0.}; //; 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 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(!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(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	}
1ff24d0a	183	if(fComputeMultiplicity) 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
54f95b99	194	Float_t lc1[3]={0.,0.,0.}; // for(Int_t ii=0; ii<3; ii++) lc1[ii]=0.;
54f95b99	195	Float_t gc1[3]={0.,0.,0.}; // ; for(Int_t ii=0; ii<3; ii++) gc1[ii]=0.;
ecc64c3f	196	// lc2 and gc2 are local and global coordinates for layer 2
54f95b99	197	Float_t lc2[3]={0.,0.,0.}; // ; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
54f95b99	198	Float_t gc2[3]={0.,0.,0.}; //; 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	*/
54f95b99	238	Int_t nEntriesMod[kNSPDMod];
	239	TClonesArray* recpArr[kNSPDMod];
	240	for(Int_t modul=0; modul<kNSPDMod; ++modul) {
	241	if(!fUseModule[modul]) {
	242	nEntriesMod[modul]=0;
	243	recpArr[modul]=0;
	244	} else {
	245	recpArr[modul]=rpcont->UncheckedGetClusters(modul);
	246	nEntriesMod[modul]=recpArr[modul]->GetEntriesFast();
	247	}
	248	}
	249
27167524	250	Int_t maxdim=TMath::Min(nrpL1*nrpL2,50000); // temporary; to limit the size in PbPb
f5f6da22	251	static TClonesArray points("AliITSZPoint",maxdim);
7299b7a8	252	Int_t nopoints =0;
27167524	253	for(Int_t modul1= fFirstL1; modul1<=fLastL1;modul1++){ // Loop on modules of layer 1
8c42830a	254	if(!fUseModule[modul1]) continue;
27167524	255	UShort_t ladder=int(modul1/4)+1; // ladders are numbered starting from 1
54f95b99	256	TClonesArray *prpl1=recpArr[modul1]; //rpcont->UncheckedGetClusters(modul1);
54f95b99	257	Int_t nrecp1 = nEntriesMod[modul1]; //prpl1->GetEntries();
27167524	258	for(Int_t j1=0;j1<nrecp1;j1++){
54f95b99	259	AliITSRecPoint recp1 = (AliITSRecPoint)prpl1->At(j1);
54f95b99	260	recp1->GetGlobalXYZ(gc1);
308c2f7c	261	gc1[0]-=GetNominalPos()[0]; // Possible beam offset in the bending plane
308c2f7c	262	gc1[1]-=GetNominalPos()[1]; // " "
27167524	263	Float_t phi1 = TMath::ATan2(gc1[1],gc1[0]);
54f95b99	264	if(phi1<0)phi1+=TMath::TwoPi();
27167524	265	for(Int_t ladl2=0 ; ladl2<fLadOnLay2*2+1;ladl2++){
	266	for(Int_t k=0;k<4;k++){
	267	Int_t ladmod=fLadders[ladder-1]+ladl2;
32e63e47	268	if(ladmod>AliITSgeomTGeo::GetNLadders(2)) ladmod=ladmod-AliITSgeomTGeo::GetNLadders(2);
32e63e47	269	Int_t modul2=AliITSgeomTGeo::GetModuleIndex(2,ladmod,k+1);
8c42830a	270	if(!fUseModule[modul2]) continue;
54f95b99	271	itsRec=recpArr[modul2]; // rpcont->UncheckedGetClusters(modul2);
54f95b99	272	Int_t nrecp2 = nEntriesMod[modul2]; // itsRec->GetEntries();
27167524	273	for(Int_t j2=0;j2<nrecp2;j2++){
54f95b99	274	AliITSRecPoint recp2 = (AliITSRecPoint)itsRec->At(j2);
54f95b99	275	recp2->GetGlobalXYZ(gc2);
308c2f7c	276	gc2[0]-=GetNominalPos()[0];
308c2f7c	277	gc2[1]-=GetNominalPos()[1];
27167524	278	Float_t phi2 = TMath::ATan2(gc2[1],gc2[0]);
54f95b99	279	if(phi2<0)phi2+=TMath::TwoPi();
27167524	280
27167524	281	Float_t diff = TMath::Abs(phi2-phi1);
54f95b99	282	if(diff>TMath::Pi())diff=TMath::TwoPi()-diff;
27167524	283	if(diff<fDiffPhiMax){
54f95b99	284	Float_t r1=TMath::Sqrt(gc1[0]gc1[0]+gc1[1]gc1[1]);
	285	Float_t zc1=gc1[2];
	286	Float_t erz1=recp1->GetSigmaZ2();
	287	Float_t r2=TMath::Sqrt(gc2[0]gc2[0]+gc2[1]gc2[1]);
	288	Float_t zc2=gc2[2];
	289	Float_t erz2=recp2->GetSigmaZ2();
27167524	290	// Float_t tgth=(zc2[j]-zc1[i])/(r2-r1); // slope (used for multiple scattering)
27167524	291	Float_t zr0=(r2zc1-r1zc2)/(r2-r1); //Z @ null radius
54f95b99	292	Float_t ezr0q=(r2r2erz1+r1r1erz2)/((r2-r1)*(r2-r1)); //error on Z @ null radius
	293	/*
	294	// Multiple scattering
	295	ezr0q+=r1r1(1+tgthtgth)theta2L1/2; // multiple scattering in layer 1
	296	ezr0q+=rBPrBP(1+tgthtgth)theta2BP/2; // multiple scattering in beam pipe
	297	*/
f5f6da22	298	if(nopoints<maxdim) new(points[nopoints++])AliITSZPoint(zr0,ezr0q);
27167524	299	fZCombc->Fill(zr0);
	300	}
	301	}
27167524	302	}
0c6af5c9	303	}
	304	}
	305	}
7299b7a8	306
f5f6da22	307	points.Sort();
7299b7a8	308
6fd990e3	309	Double_t contents = fZCombc->GetEntries()- fZCombc->GetBinContent(0)-fZCombc->GetBinContent(nbincoarse+1);
	310	if(contents<1.){
	311	// Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
0c6af5c9	312	ResetHistograms();
6fd990e3	313	fCurrentVertex = new AliESDVertex(0.,5.3,-1);
f5f6da22	314	points.Clear();
6fd990e3	315	return;
0c6af5c9	316	}
ecc64c3f	317
ecc64c3f	318	TH1F *hc = fZCombc;
ecc64c3f	319
ecc64c3f	320
8c42830a	321	if(hc->GetBinContent(hc->GetMaximumBin())<3)hc->Rebin(4);
debb3dde	322	Int_t binmin,binmax;
debb3dde	323	Int_t nPeaks=GetPeakRegion(hc,binmin,binmax);
4be14da5	324	if(nPeaks==2)AliDebug(2,"2 peaks found");
7299b7a8	325	Float_t zm =0.;
	326	Float_t ezm =0.;
	327	Float_t lim1 = hc->GetBinLowEdge(binmin);
	328	Float_t lim2 = hc->GetBinLowEdge(binmax)+hc->GetBinWidth(binmax);
4be14da5	329	Float_t widthSR=lim2-lim1;
7299b7a8	330
	331	if(nPeaks ==1 && (lim2-lim1)<fWindowWidth){
	332	Float_t c=(lim1+lim2)/2.;
	333	lim1=c-fWindowWidth/2.;
	334	lim2=c+fWindowWidth/2.;
ecc64c3f	335	}
7299b7a8	336	Int_t niter = 0, ncontr=0;
	337	do {
	338	// symmetrization
	339	if(zm !=0.){
	340	Float_t semilarg=TMath::Min((lim2-zm),(zm-lim1));
	341	lim1=zm - semilarg;
	342	lim2=zm + semilarg;
0c6af5c9	343	}
7299b7a8	344
	345	zm=0.;
	346	ezm=0.;
	347	ncontr=0;
f5f6da22	348	for(Int_t i =0; i<points.GetEntries(); i++){
f5f6da22	349	AliITSZPoint* p=(AliITSZPoint*)points.UncheckedAt(i);
7299b7a8	350	if(p->GetZ()>lim1 && p->GetZ()<lim2){
	351	Float_t deno = p->GetErrZ();
	352	zm+=p->GetZ()/deno;
	353	ezm+=1./deno;
	354	ncontr++;
ecc64c3f	355	}
ecc64c3f	356	}
2b8b4286	357	if(ezm>0) {
	358	zm/=ezm;
	359	ezm=TMath::Sqrt(1./ezm);
	360	}
7299b7a8	361	niter++;
7299b7a8	362	} while(niter<10 && TMath::Abs((zm-lim1)-(lim2-zm))>fTolerance);
4be14da5	363	if(nPeaks==2) ezm=widthSR;
7299b7a8	364	fCurrentVertex = new AliESDVertex(zm,ezm,ncontr);
12e3ead8	365	fCurrentVertex->SetTitle("vertexer: Z");
bf2e0ad4	366	fCurrentVertex->SetDispersion(fDiffPhiMax);
6b4d9537	367	fNoVertices=1;
f5f6da22	368	points.Clear();
8c32ba44	369	if(ncontr>fMinTrackletsForPilup){
	370	Float_t secPeakPos;
	371	Int_t ncontr2=FindSecondPeak(fZCombc,binmin,binmax,secPeakPos);
	372	if(ncontr2>=fMinTrackletsForPilup){
	373	fIsPileup=kTRUE;
6b4d9537	374	fNoVertices=2;
8c32ba44	375	fZpuv=secPeakPos;
8c32ba44	376	fNTrpuv=ncontr2;
6b4d9537	377	AliESDVertex secondVert(secPeakPos,0.1,ncontr2);
	378	fVertArray = new AliESDVertex[2];
	379	fVertArray[0]=(*fCurrentVertex);
	380	fVertArray[1]=secondVert;
8c32ba44	381	}
2c3a7cd0	382	}
6b4d9537	383	if(fNoVertices==1){
	384	fVertArray = new AliESDVertex[1];
	385	fVertArray[0]=(*fCurrentVertex);
	386	}
	387
0c6af5c9	388	ResetHistograms();
6fd990e3	389	return;
0c6af5c9	390	}
0c6af5c9	391
8c32ba44	392	//_____________________________________________________________________
	393	Int_t AliITSVertexerZ::FindSecondPeak(TH1F* h, Int_t binmin,Int_t binmax, Float_t& secPeakPos){
	394	for(Int_t i=binmin-1;i<=binmax+1;i++){
	395	h->SetBinContent(i,0.);
	396	}
	397	Int_t secPeakBin=h->GetMaximumBin();
	398	secPeakPos=h->GetBinCenter(secPeakBin);
70a3283c	399	Int_t secPeakCont=(Int_t)h->GetBinContent(secPeakBin);
	400	secPeakCont+=(Int_t)h->GetBinContent(secPeakBin-1);
	401	secPeakCont+=(Int_t)h->GetBinContent(secPeakBin+1);
	402	secPeakCont+=(Int_t)h->GetBinContent(secPeakBin-2);
	403	secPeakCont+=(Int_t)h->GetBinContent(secPeakBin+2);
8c32ba44	404	return secPeakCont;
	405	}
	406
0c6af5c9	407	//_____________________________________________________________________
	408	void AliITSVertexerZ::ResetHistograms(){
	409	// delete TH1 data members
	410	if(fZCombc)delete fZCombc;
0c6af5c9	411	fZCombc = 0;
0c6af5c9	412	}
0c6af5c9	413
0c6af5c9	414	//________________________________________________________
	415	void AliITSVertexerZ::PrintStatus() const {
	416	// Print current status
	417	cout <<"=======================================================\n";
	418	cout <<" First layer first and last modules: "<<fFirstL1<<", ";
	419	cout <<fLastL1<<endl;
	420	cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
	421	cout <<fLastL2<<endl;
	422	cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
	423	cout <<"Limits for Z histograms: "<<fLowLim<<"; "<<fHighLim<<endl;
7299b7a8	424	cout <<"Bin sizes for coarse z histos "<<fStepCoarse<<endl;
0c6af5c9	425	cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
ecc64c3f	426	if(fZCombc){
	427	cout<<"fZCombc exists - entries="<<fZCombc->GetEntries()<<endl;
	428	}
	429	else{
	430	cout<<"fZCombc does not exist\n";
	431	}
0c6af5c9	432
	433	cout <<"=======================================================\n";
	434	}
	435