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