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