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