[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.05,0.1,0.2,0.3,0.5};
  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;

  // By default fFirstL1=0 and fLastL1=79
  for(Int_t module= fFirstL1; module<=fLastL1;module++){
    itsRec=rpcont->UncheckedGetClusters(module);
    nrpL1+= itsRec->GetEntries();
  }
  //By default fFirstL2=80 and fLastL2=239
  for(Int_t module= fFirstL2; module<=fLastL2;module++){
    itsRec=rpcont->UncheckedGetClusters(module);
    nrpL2+= itsRec->GetEntries();
  }
  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)AliWarning("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);

  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);
  fCurrentVertex = new AliESDVertex(zm,ezm,ncontr);
  fCurrentVertex->SetTitle("vertexer: Z");
  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
	109	Float_t defaults[5]={0.05,0.1,0.2,0.3,0.5};
	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;
27167524	210
ecc64c3f	211	// By default fFirstL1=0 and fLastL1=79
0c6af5c9	212	for(Int_t module= fFirstL1; module<=fLastL1;module++){
b21c1af0	213	itsRec=rpcont->UncheckedGetClusters(module);
0c6af5c9	214	nrpL1+= itsRec->GetEntries();
0c6af5c9	215	}
ecc64c3f	216	//By default fFirstL2=80 and fLastL2=239
0c6af5c9	217	for(Int_t module= fFirstL2; module<=fLastL2;module++){
b21c1af0	218	itsRec=rpcont->UncheckedGetClusters(module);
0c6af5c9	219	nrpL2+= itsRec->GetEntries();
0c6af5c9	220	}
0c6af5c9	221	if(nrpL1 == 0 \|\| nrpL2 == 0){
fea0ad94	222	AliDebug(1,Form("No RecPoints in at least one SPD layer (%d %d)",nrpL1,nrpL2));
0c6af5c9	223	ResetHistograms();
6fd990e3	224	fCurrentVertex = new AliESDVertex(0.,5.3,-2);
6fd990e3	225	return;
0c6af5c9	226	}
ecc64c3f	227	// Force a coarse bin size of 200 microns if the number of clusters on layer 2
	228	// is low
	229	if(nrpL2<fPPsetting[0])SetBinWidthCoarse(fPPsetting[1]);
ecc64c3f	230	// By default nbincoarse=(10+10)/0.01=2000
ecc64c3f	231	Int_t nbincoarse = static_cast<Int_t>((fHighLim-fLowLim)/fStepCoarse);
ecc64c3f	232	if(fZCombc)delete fZCombc;
ecc64c3f	233	fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);
7299b7a8	234
27167524	235	/* Test the ffect of mutiple scatternig on error. Negligible
7299b7a8	236	// Multiple scattering
	237	Float_t beta=1.,pmed=0.875; //pmed=875 MeV (for tracks with dphi<0.01 rad)
	238	Float_t beta2=beta*beta;
	239	Float_t p2=pmed*pmed;
	240	Float_t rBP=3; //Beam Pipe radius = 3cm
	241	Float_t dBP=0.08/35.3; // 800 um of Be
	242	Float_t dL1=0.01; //approx. 1% of radiation length
	243	Float_t theta2BP=14.114.1/(beta2p21e6)TMath::Abs(dBP);
	244	Float_t theta2L1=14.114.1/(beta2p21e6)TMath::Abs(dL1);
	245	*/
27167524	246	Int_t maxdim=TMath::Min(nrpL1*nrpL2,50000); // temporary; to limit the size in PbPb
f5f6da22	247	static TClonesArray points("AliITSZPoint",maxdim);
7299b7a8	248	Int_t nopoints =0;
27167524	249	for(Int_t modul1= fFirstL1; modul1<=fLastL1;modul1++){ // Loop on modules of layer 1
8c42830a	250	if(!fUseModule[modul1]) continue;
27167524	251	UShort_t ladder=int(modul1/4)+1; // ladders are numbered starting from 1
b21c1af0	252	TClonesArray *prpl1=rpcont->UncheckedGetClusters(modul1);
b21c1af0	253	Int_t nrecp1 = prpl1->GetEntries();
27167524	254	for(Int_t j1=0;j1<nrecp1;j1++){
b21c1af0	255	AliITSRecPoint recp = (AliITSRecPoint)prpl1->At(j1);
32e63e47	256	recp->GetGlobalXYZ(gc1);
308c2f7c	257	gc1[0]-=GetNominalPos()[0]; // Possible beam offset in the bending plane
308c2f7c	258	gc1[1]-=GetNominalPos()[1]; // " "
27167524	259	Float_t r1=TMath::Sqrt(gc1[0]gc1[0]+gc1[1]gc1[1]);
	260	Float_t phi1 = TMath::ATan2(gc1[1],gc1[0]);
	261	if(phi1<0)phi1+=2*TMath::Pi();
	262	Float_t zc1=gc1[2];
	263	Float_t erz1=recp->GetSigmaZ2();
	264	for(Int_t ladl2=0 ; ladl2<fLadOnLay2*2+1;ladl2++){
	265	for(Int_t k=0;k<4;k++){
	266	Int_t ladmod=fLadders[ladder-1]+ladl2;
32e63e47	267	if(ladmod>AliITSgeomTGeo::GetNLadders(2)) ladmod=ladmod-AliITSgeomTGeo::GetNLadders(2);
32e63e47	268	Int_t modul2=AliITSgeomTGeo::GetModuleIndex(2,ladmod,k+1);
8c42830a	269	if(!fUseModule[modul2]) continue;
b21c1af0	270	itsRec=rpcont->UncheckedGetClusters(modul2);
27167524	271	Int_t nrecp2 = itsRec->GetEntries();
	272	for(Int_t j2=0;j2<nrecp2;j2++){
	273	recp = (AliITSRecPoint*)itsRec->At(j2);
32e63e47	274	recp->GetGlobalXYZ(gc2);
308c2f7c	275	gc2[0]-=GetNominalPos()[0];
308c2f7c	276	gc2[1]-=GetNominalPos()[1];
27167524	277	Float_t r2=TMath::Sqrt(gc2[0]gc2[0]+gc2[1]gc2[1]);
	278	Float_t phi2 = TMath::ATan2(gc2[1],gc2[0]);
	279	if(phi2<0)phi2+=2*TMath::Pi();
	280	Float_t zc2=gc2[2];
	281	Float_t erz2=recp->GetSigmaZ2();
	282
	283	Float_t diff = TMath::Abs(phi2-phi1);
	284	if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
	285	if(diff<fDiffPhiMax){
	286	// Float_t tgth=(zc2[j]-zc1[i])/(r2-r1); // slope (used for multiple scattering)
	287	Float_t zr0=(r2zc1-r1zc2)/(r2-r1); //Z @ null radius
	288	Float_t ezr0q=(r2r2erz1+r1r1erz2)/(r2-r1)/(r2-r1); //error on Z @ null radius
	289	/*
	290	// Multiple scattering
7299b7a8	291	ezr0q+=r1r1(1+tgthtgth)theta2L1/2; // multiple scattering in layer 1
	292	ezr0q+=rBPrBP(1+tgthtgth)theta2BP/2; // multiple scattering in beam pipe
	293	*/
f5f6da22	294	if(nopoints<maxdim) new(points[nopoints++])AliITSZPoint(zr0,ezr0q);
27167524	295	fZCombc->Fill(zr0);
	296	}
	297	}
27167524	298	}
0c6af5c9	299	}
	300	}
	301	}
7299b7a8	302
f5f6da22	303	points.Sort();
7299b7a8	304
6fd990e3	305	Double_t contents = fZCombc->GetEntries()- fZCombc->GetBinContent(0)-fZCombc->GetBinContent(nbincoarse+1);
	306	if(contents<1.){
	307	// Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
0c6af5c9	308	ResetHistograms();
6fd990e3	309	fCurrentVertex = new AliESDVertex(0.,5.3,-1);
f5f6da22	310	points.Clear();
6fd990e3	311	return;
0c6af5c9	312	}
ecc64c3f	313
ecc64c3f	314	TH1F *hc = fZCombc;
ecc64c3f	315
ecc64c3f	316
8c42830a	317	if(hc->GetBinContent(hc->GetMaximumBin())<3)hc->Rebin(4);
debb3dde	318	Int_t binmin,binmax;
	319	Int_t nPeaks=GetPeakRegion(hc,binmin,binmax);
	320	if(nPeaks==2)AliWarning("2 peaks found");
7299b7a8	321	Float_t zm =0.;
	322	Float_t ezm =0.;
	323	Float_t lim1 = hc->GetBinLowEdge(binmin);
	324	Float_t lim2 = hc->GetBinLowEdge(binmax)+hc->GetBinWidth(binmax);
	325
	326	if(nPeaks ==1 && (lim2-lim1)<fWindowWidth){
	327	Float_t c=(lim1+lim2)/2.;
	328	lim1=c-fWindowWidth/2.;
	329	lim2=c+fWindowWidth/2.;
ecc64c3f	330	}
7299b7a8	331	Int_t niter = 0, ncontr=0;
	332	do {
	333	// symmetrization
	334	if(zm !=0.){
	335	Float_t semilarg=TMath::Min((lim2-zm),(zm-lim1));
	336	lim1=zm - semilarg;
	337	lim2=zm + semilarg;
0c6af5c9	338	}
7299b7a8	339
	340	zm=0.;
	341	ezm=0.;
	342	ncontr=0;
f5f6da22	343	for(Int_t i =0; i<points.GetEntries(); i++){
f5f6da22	344	AliITSZPoint* p=(AliITSZPoint*)points.UncheckedAt(i);
7299b7a8	345	if(p->GetZ()>lim1 && p->GetZ()<lim2){
	346	Float_t deno = p->GetErrZ();
	347	zm+=p->GetZ()/deno;
	348	ezm+=1./deno;
	349	ncontr++;
ecc64c3f	350	}
ecc64c3f	351	}
2b8b4286	352	if(ezm>0) {
	353	zm/=ezm;
	354	ezm=TMath::Sqrt(1./ezm);
	355	}
7299b7a8	356	niter++;
	357	} while(niter<10 && TMath::Abs((zm-lim1)-(lim2-zm))>fTolerance);
	358	fCurrentVertex = new AliESDVertex(zm,ezm,ncontr);
12e3ead8	359	fCurrentVertex->SetTitle("vertexer: Z");
6b4d9537	360	fNoVertices=1;
f5f6da22	361	points.Clear();
8c32ba44	362	if(ncontr>fMinTrackletsForPilup){
	363	Float_t secPeakPos;
	364	Int_t ncontr2=FindSecondPeak(fZCombc,binmin,binmax,secPeakPos);
	365	if(ncontr2>=fMinTrackletsForPilup){
	366	fIsPileup=kTRUE;
6b4d9537	367	fNoVertices=2;
8c32ba44	368	fZpuv=secPeakPos;
8c32ba44	369	fNTrpuv=ncontr2;
6b4d9537	370	AliESDVertex secondVert(secPeakPos,0.1,ncontr2);
	371	fVertArray = new AliESDVertex[2];
	372	fVertArray[0]=(*fCurrentVertex);
	373	fVertArray[1]=secondVert;
8c32ba44	374	}
2c3a7cd0	375	}
6b4d9537	376	if(fNoVertices==1){
	377	fVertArray = new AliESDVertex[1];
	378	fVertArray[0]=(*fCurrentVertex);
	379	}
	380
0c6af5c9	381	ResetHistograms();
6fd990e3	382	return;
0c6af5c9	383	}
0c6af5c9	384
8c32ba44	385	//_____________________________________________________________________
	386	Int_t AliITSVertexerZ::FindSecondPeak(TH1F* h, Int_t binmin,Int_t binmax, Float_t& secPeakPos){
	387	for(Int_t i=binmin-1;i<=binmax+1;i++){
	388	h->SetBinContent(i,0.);
	389	}
	390	Int_t secPeakBin=h->GetMaximumBin();
	391	secPeakPos=h->GetBinCenter(secPeakBin);
70a3283c	392	Int_t secPeakCont=(Int_t)h->GetBinContent(secPeakBin);
	393	secPeakCont+=(Int_t)h->GetBinContent(secPeakBin-1);
	394	secPeakCont+=(Int_t)h->GetBinContent(secPeakBin+1);
	395	secPeakCont+=(Int_t)h->GetBinContent(secPeakBin-2);
	396	secPeakCont+=(Int_t)h->GetBinContent(secPeakBin+2);
8c32ba44	397	return secPeakCont;
	398	}
	399
0c6af5c9	400	//_____________________________________________________________________
	401	void AliITSVertexerZ::ResetHistograms(){
	402	// delete TH1 data members
	403	if(fZCombc)delete fZCombc;
0c6af5c9	404	fZCombc = 0;
0c6af5c9	405	}
0c6af5c9	406
0c6af5c9	407	//________________________________________________________
	408	void AliITSVertexerZ::PrintStatus() const {
	409	// Print current status
	410	cout <<"=======================================================\n";
	411	cout <<" First layer first and last modules: "<<fFirstL1<<", ";
	412	cout <<fLastL1<<endl;
	413	cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
	414	cout <<fLastL2<<endl;
	415	cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
	416	cout <<"Limits for Z histograms: "<<fLowLim<<"; "<<fHighLim<<endl;
7299b7a8	417	cout <<"Bin sizes for coarse z histos "<<fStepCoarse<<endl;
0c6af5c9	418	cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
ecc64c3f	419	if(fZCombc){
	420	cout<<"fZCombc exists - entries="<<fZCombc->GetEntries()<<endl;
	421	}
	422	else{
	423	cout<<"fZCombc does not exist\n";
	424	}
0c6af5c9	425
	426	cout <<"=======================================================\n";
	427	}
	428