[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<TH1.h>
#include <TString.h>
#include<TTree.h>
#include "AliITSLoader.h"
#include "AliITSgeom.h"
#include "AliITSDetTypeRec.h"
#include "AliITSRecPoint.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),
fX0(0.),
fY0(0.),
fZFound(0),
fZsig(0.),
fZCombc(0),
fZCombv(0),
fZCombf(0),
fLowLim(0.),
fHighLim(0.),
fStepCoarse(0),
fStepFine(0),
fTolerance(0.),
fMaxIter(0){
  // Default constructor
  SetDiffPhiMax(0);
  SetFirstLayerModules(0);
  SetSecondLayerModules(0);
  SetLowLimit(0.);
  SetHighLimit(0.);
  SetBinWidthCoarse(0.);
  SetBinWidthFine(0.);
  SetTolerance(0.);
  SetPPsetting(0.,0.);
  ConfigIterations();
}

//______________________________________________________________________
AliITSVertexerZ::AliITSVertexerZ(TString fn, Float_t x0, Float_t y0):AliITSVertexer(fn),
fFirstL1(0),
fLastL1(0),
fFirstL2(0),
fLastL2(0),
fDiffPhiMax(0),
fX0(x0),
fY0(y0),
fZFound(0),
fZsig(0.),
fZCombc(0),
fZCombv(0),
fZCombf(0),
fLowLim(0.),
fHighLim(0.),
fStepCoarse(0),
fStepFine(0),
fTolerance(0.),
fMaxIter(0) {
  // Standard Constructor
  SetDiffPhiMax();
  SetFirstLayerModules();
  SetSecondLayerModules();
  SetLowLimit();
  SetHighLimit();
  SetBinWidthCoarse();
  SetBinWidthFine();
  SetTolerance();
  SetPPsetting();
  ConfigIterations();

}

//______________________________________________________________________
AliITSVertexerZ::AliITSVertexerZ(const AliITSVertexerZ &vtxr) : AliITSVertexer(vtxr),
fFirstL1(vtxr.fFirstL1),
fLastL1(vtxr.fLastL1),
fFirstL2(vtxr.fFirstL2),
fLastL2(vtxr.fLastL2),
fDiffPhiMax(vtxr.fDiffPhiMax),
fX0(vtxr.fX0),
fY0(vtxr.fY0),
fZFound(vtxr.fZFound),
fZsig(vtxr.fZsig),
fZCombc(vtxr.fZCombc),
fZCombv(vtxr.fZCombv),
fZCombf(vtxr.fZCombf),
fLowLim(vtxr.fLowLim),
fHighLim(vtxr.fHighLim),
fStepCoarse(vtxr.fStepCoarse),
fStepFine(vtxr.fStepFine),
fTolerance(vtxr.fTolerance),
fMaxIter(vtxr.fMaxIter) {
  // Copy constructor

}

//______________________________________________________________________
AliITSVertexerZ& AliITSVertexerZ::operator=(const AliITSVertexerZ&  vtxr ){
  // Assignment operator

  this->~AliITSVertexerZ();
  new(this) AliITSVertexerZ(vtxr);
  return *this;
}

//______________________________________________________________________
AliITSVertexerZ::~AliITSVertexerZ() {
  // Destructor
  delete fZCombc;
  delete fZCombf;
  delete fZCombv;
}

//______________________________________________________________________
void AliITSVertexerZ::ConfigIterations(Int_t noiter,Float_t *ptr){
  // configure the iterative procedure to gain efficiency for
  // pp events with very low multiplicity
  Float_t defaults[5]={0.05,0.1,0.2,0.3,0.5};
  fMaxIter=noiter;
  if(noiter>5){
    Error("ConfigIterations","Maximum number of iterations is 5\n");
    fMaxIter=5;
  }
  for(Int_t j=0;j<5;j++)fPhiDiffIter[j]=defaults[j];
  if(ptr)for(Int_t j=0;j<fMaxIter;j++)fPhiDiffIter[j]=ptr[j];
}

//______________________________________________________________________
Int_t AliITSVertexerZ::GetPeakRegion(TH1F*h, Int_t &binmin, Int_t &binmax) const {
  // 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(Int_t evnumber){
  // Defines the AliESDVertex for the current event
  VertexZFinder(evnumber);
  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(evnumber);
      fDiffPhiMax=diffPhiMaxOrig;
    }
  }
  FindMultiplicity(evnumber);
  return fCurrentVertex;
}  


//______________________________________________________________________
void AliITSVertexerZ::VertexZFinder(Int_t evnumber){
  // Defines the AliESDVertex for the current event
  fCurrentVertex = 0;
  AliRunLoader *rl =AliRunLoader::GetRunLoader();
  AliITSLoader* itsLoader = (AliITSLoader*)rl->GetLoader("ITSLoader");
  AliITSgeom* geom = itsLoader->GetITSgeom();
  itsLoader->LoadRecPoints();
  rl->GetEvent(evnumber);

  AliITSDetTypeRec detTypeRec;

  TTree *tR = itsLoader->TreeR();
  detTypeRec.SetTreeAddressR(tR);
  TClonesArray *itsRec  = 0;
  // lc and gc are local and global coordinates for layer 1
  Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
  Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
  // lc2 and gc2 are local and global coordinates for layer 2
  Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
  Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;

  itsRec = detTypeRec.RecPoints();
  TBranch *branch;
  branch = tR->GetBranch("ITSRecPoints");

  Int_t nrpL1 = 0;
  Int_t nrpL2 = 0;
  // By default fFirstL1=0 and fLastL1=79
  // This loop counts the number of recpoints on layer1 (central modules)
  for(Int_t module= fFirstL1; module<=fLastL1;module++){
    // Keep only central modules
    if(module%4==0 || module%4==3)continue;
    //   cout<<"Procesing module "<<module<<" ";
    branch->GetEvent(module);
    //    cout<<"Number of clusters "<<clusters->GetEntries()<<endl;
    nrpL1+= itsRec->GetEntries();
    detTypeRec.ResetRecPoints();
  }
  //By default fFirstL2=80 and fLastL2=239
  //This loop counts the number of RP on layer 2
  for(Int_t module= fFirstL2; module<=fLastL2;module++){
    branch->GetEvent(module);
    nrpL2+= itsRec->GetEntries();
    detTypeRec.ResetRecPoints();
  }
  if(nrpL1 == 0 || nrpL2 == 0){
    ResetHistograms();
    itsLoader->UnloadRecPoints();
    fCurrentVertex = new AliESDVertex(0.,5.3,-2);
    return;
  }
  // The vertex finding is attempted only if the number of RP is !=0 on
  // both layers
  Float_t *xc1 = new Float_t [nrpL1]; // coordinates of the L1 Recpoints
  Float_t *yc1 = new Float_t [nrpL1];
  Float_t *zc1 = new Float_t [nrpL1];
  Float_t *phi1 = new Float_t [nrpL1];
  Float_t *xc2 = new Float_t [nrpL2]; // coordinates of the L1 Recpoints
  Float_t *yc2 = new Float_t [nrpL2];
  Float_t *zc2 = new Float_t [nrpL2];
  Float_t *phi2 = new Float_t [nrpL2];
  Int_t ind = 0;// running index for RP
  // 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 nbinfine = (10+10)/0.0005=40000
  Int_t nbinfine = static_cast<Int_t>((fHighLim-fLowLim)/fStepFine);
  // By default nbincoarse=(10+10)/0.01=2000
  Int_t nbincoarse = static_cast<Int_t>((fHighLim-fLowLim)/fStepCoarse);
  // Set stepverycoarse = 3*fStepCoarse
  Int_t nbinvcoarse = static_cast<Int_t>((fHighLim-fLowLim)/(3.*fStepCoarse));
  if(fZCombc)delete fZCombc;
  fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);
  if(fZCombv)delete fZCombv;
  fZCombv = new TH1F("fZCombv","Z",nbinvcoarse,fLowLim,fLowLim+nbinvcoarse*3.*fStepCoarse);
  if(fZCombf)delete fZCombf;
  fZCombf = new TH1F("fZCombf","Z",nbinfine,fLowLim,fLowLim+nbinfine*fStepFine);
  // Loop on modules of layer 1 (restricted to central modules)
  for(Int_t module= fFirstL1; module<=fLastL1;module++){
    if(module%4==0 || module%4==3)continue;
    branch->GetEvent(module);
    Int_t nrecp1 = itsRec->GetEntries();
    for(Int_t j=0;j<nrecp1;j++){
      AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
      // Local coordinates of this recpoint
      lc[0]=recp->GetDetLocalX();
      lc[2]=recp->GetDetLocalZ();
      geom->LtoG(module,lc,gc);
      // Global coordinates of this recpoints
      gc[0]-=fX0; // Possible beam offset in the bending plane
      gc[1]-=fY0; //   "               "
      xc1[ind]=gc[0];
      yc1[ind]=gc[1];
      zc1[ind]=gc[2];
      // azimuthal angle is computed in the interval 0 --> 2*pi
      phi1[ind] = TMath::ATan2(gc[1],gc[0]);
      if(phi1[ind]<0)phi1[ind]=2*TMath::Pi()+phi1[ind];
      ind++;
    }
    detTypeRec.ResetRecPoints();
  }
  ind = 0; // the running index is reset for Layer 2
  for(Int_t module= fFirstL2; module<=fLastL2;module++){
    branch->GetEvent(module);
    Int_t nrecp2 = itsRec->GetEntries();
    for(Int_t j=0;j<nrecp2;j++){
      AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
      lc[0]=recp->GetDetLocalX();
      lc[2]=recp->GetDetLocalZ();
      geom->LtoG(module,lc,gc);
      gc[0]-=fX0;
      gc[1]-=fY0;
      xc2[ind]=gc[0];
      yc2[ind]=gc[1];
      zc2[ind]=gc[2];
      phi2[ind] = TMath::ATan2(gc[1],gc[0]);
      if(phi2[ind]<0)phi2[ind]=2*TMath::Pi()+phi2[ind];
      ind++;
    }
    detTypeRec.ResetRecPoints();
  }
 
  // Int_t nolines=0;
  for(Int_t i=0;i<nrpL1;i++){ // loop on L1 RP
    Float_t r1=TMath::Sqrt(xc1[i]*xc1[i]+yc1[i]*yc1[i]); // radius L1 RP
    for(Int_t j=0;j<nrpL2;j++){ // loop on L2 RP
      Float_t diff = TMath::Abs(phi2[j]-phi1[i]); // diff in azimuth
      if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff; //diff<pi
      if(diff<fDiffPhiMax){ // cut on 10 milliradians by def.
	Float_t r2=TMath::Sqrt(xc2[j]*xc2[j]+yc2[j]*yc2[j]); // radius L2 RP
	Float_t zr0=(r2*zc1[i]-r1*zc2[j])/(r2-r1); //Z @ null radius
	fZCombf->Fill(zr0);
	fZCombc->Fill(zr0);
	fZCombv->Fill(zr0);
	Double_t pA[3];
	Double_t pB[3];
	pA[0]=xc1[i];
	pA[1]=yc1[i];
	pA[2]=zc1[i];
	pB[0]=xc2[j];
	pB[1]=yc2[j];
	pB[2]=zc2[j];
	//	MakeTracklet(pA,pB,nolines);
      }
    }
  }
  delete [] xc1;
  delete [] yc1;
  delete [] zc1;
  delete [] phi1;
  delete [] xc2;
  delete [] yc2;
  delete [] zc2;
  delete [] phi2;
  Double_t contents = fZCombc->GetEntries()- fZCombc->GetBinContent(0)-fZCombc->GetBinContent(nbincoarse+1);
  if(contents<1.){
    //    Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
    ResetHistograms();
    itsLoader->UnloadRecPoints();
    fCurrentVertex = new AliESDVertex(0.,5.3,-1);
    return;
  }

  TH1F *hc = fZCombc;

  
  if(hc->GetBinContent(hc->GetMaximumBin())<3)hc = fZCombv;
  Int_t binmin,binmax;
  Int_t nPeaks=GetPeakRegion(hc,binmin,binmax);   
  if(nPeaks==2)AliWarning("2 peaks found");
  Int_t ii1=1+static_cast<Int_t>((hc->GetBinLowEdge(binmin)-fLowLim)/fStepFine);
  Int_t ii2=1+static_cast<Int_t>((hc->GetBinLowEdge(binmax)+hc->GetBinWidth(binmax)-fLowLim)/fStepFine);
  Float_t centre = 0.;
  Int_t nn=0;
  for(Int_t ii=ii1;ii<=ii2;ii++){
    centre+= fZCombf->GetBinCenter(ii)*fZCombf->GetBinContent(ii);
    nn+=static_cast<Int_t>(fZCombf->GetBinContent(ii));
  }
  if (nn) centre/=nn;

  // n1 is the bin number of fine histogram containing the point located 1 coarse bin less than "centre"
  Int_t n1 = 1+static_cast<Int_t>((centre-hc->GetBinWidth(1)-fLowLim)/fStepFine);
  // n2 is the bin number of fine histogram containing the point located 1 coarse bin more than "centre"
  Int_t n2 = 1+static_cast<Int_t>((centre+hc->GetBinWidth(1)-fLowLim)/fStepFine);
  if(n1<1)n1=1;
  if(n2>nbinfine)n2=nbinfine;
  Int_t niter = 0;
  Bool_t goon = kTRUE;
  Int_t num=0;
  Bool_t last = kFALSE;

  while(goon || last){
    fZFound = 0.;
    fZsig = 0.;
    num=0;
    // at the end of the loop:
    // fZFound = N*(Average Z) - where N is the number of tracklets
    // num=N
    // fZsig = N*Q - where Q is the average of Z*Z
    for(Int_t n=n1;n<=n2;n++){
      fZFound+=fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
      num+=static_cast<Int_t>(fZCombf->GetBinContent(n));
      fZsig+=fZCombf->GetBinCenter(n)*fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
    }
    if(num<2){
      fZsig = 5.3;  // Default error from the beam sigmoid
    }
    else{
      Float_t radi =  fZsig/(num-1)-fZFound*fZFound/num/(num-1);
      // radi = square root of sample variance of Z
      if(radi>0.)fZsig=TMath::Sqrt(radi);
      else fZsig=5.3;  // Default error from the beam sigmoid
      // fZfound - Average Z
      fZFound/=num;
    }
    if(!last){
      // goon is true if the distance between the found Z and the lower bin differs from the distance between the found Z and
      // the upper bin by more than tolerance (0.002)
      goon = TMath::Abs(TMath::Abs(fZFound-fZCombf->GetBinCenter(n1))-TMath::Abs(fZFound-fZCombf->GetBinCenter(n2)))>fTolerance;
      // a window in the fine grained histogram is centered aroung the found Z. The width is 2 bins of
      // the coarse grained histogram
      if(num>0){
	n1 = 1+static_cast<Int_t>((fZFound-hc->GetBinWidth(1)-fLowLim)/fStepFine);
	if(n1<1)n1=1;
	n2 = 1+static_cast<Int_t>((fZFound+hc->GetBinWidth(1)-fLowLim)/fStepFine);
	if(n2>nbinfine)n2=nbinfine;
      }
      else {
	n1 = 1+static_cast<Int_t>((centre-(niter+2)*hc->GetBinWidth(1)-fLowLim)/fStepFine);
	n2 = 1+static_cast<Int_t>((centre+(niter+2)*hc->GetBinWidth(1)-fLowLim)/fStepFine);
	if(n1<1)n1=1;
	if(n2>nbinfine)n2=nbinfine;
      }
      niter++;
      // no more than 10 adjusting iterations
      if(niter>=10)goon = kFALSE;

      if((fZsig> 0.0001) && !goon && num>5){
	last = kTRUE;
	n1 = 1+static_cast<Int_t>((fZFound-fZsig-fLowLim)/fStepFine);
	if(n1<1)n1=1;
	n2 = 1+static_cast<Int_t>((fZFound+fZsig-fLowLim)/fStepFine);  
	if(n2>nbinfine)n2=nbinfine;
      }
    }
    else {
      last = kFALSE;
    }

  }
  //  if(fDebug>0)cout<<"Numer of Iterations "<<niter<<endl<<endl;
  //  if(num!=0)fZsig/=TMath::Sqrt(num);
  if (fZsig<=0) fZsig=5.3; // Default error from the beam sigmoid
  fCurrentVertex = new AliESDVertex(fZFound,fZsig,num);
  fCurrentVertex->SetTitle("vertexer: B");
  ResetHistograms();
  itsLoader->UnloadRecPoints();
  return;
}

//_____________________________________________________________________
void AliITSVertexerZ::ResetHistograms(){
  // delete TH1 data members
  if(fZCombc)delete fZCombc;
  if(fZCombf)delete fZCombf;
  if(fZCombv)delete fZCombv;
  fZCombc = 0;
  fZCombf = 0;
  fZCombv = 0;
}

//______________________________________________________________________
void AliITSVertexerZ::FindVertices(){
  // computes the vertices of the events in the range FirstEvent - LastEvent
  AliRunLoader *rl = AliRunLoader::GetRunLoader();
  AliITSLoader* itsLoader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
  itsLoader->ReloadRecPoints();
  for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
    //  cout<<"Processing event "<<i<<endl;
    rl->GetEvent(i);
    FindVertexForCurrentEvent(i);
    if(fCurrentVertex){
      WriteCurrentVertex();
    }
    /*
    else {
      if(fDebug>0){
	cout<<"Vertex not found for event "<<i<<endl;
	cout<<"fZFound = "<<fZFound<<", fZsig= "<<fZsig<<endl;
      }
    }
    */
  }
}

//________________________________________________________
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 and fine z histos "<<fStepCoarse<<"; "<<fStepFine<<endl;
  cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
  cout <<" Debug flag: "<<fDebug<<endl;
  cout <<"First event to be processed "<<fFirstEvent;
  cout <<"\n Last event to be processed "<<fLastEvent<<endl;
  if(fZCombc){
    cout<<"fZCombc exists - entries="<<fZCombc->GetEntries()<<endl;
  }
  else{
    cout<<"fZCombc does not exist\n";
  }
  if(fZCombv){
    cout<<"fZCombv exists - entries="<<fZCombv->GetEntries()<<endl;
  }
  else{
    cout<<"fZCombv does not exist\n";
  }
   if(fZCombf){
    cout<<"fZCombf exists - entries="<<fZCombv->GetEntries()<<endl;
  }
  else{
    cout<<"fZCombf 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>
0c6af5c9	17	#include<TH1.h>
d531b7a1	18	#include <TString.h>
0c6af5c9	19	#include<TTree.h>
d531b7a1	20	#include "AliITSLoader.h"
d531b7a1	21	#include "AliITSgeom.h"
7d62fb64	22	#include "AliITSDetTypeRec.h"
d531b7a1	23	#include "AliITSRecPoint.h"
0c6af5c9	24
	25	/////////////////////////////////////////////////////////////////
	26	// this class implements a fast method to determine
	27	// the Z coordinate of the primary vertex
	28	// for p-p collisions it seems to give comparable or better results
	29	// with respect to what obtained with AliITSVertexerPPZ
	30	// It can be used successfully with Pb-Pb collisions
	31	////////////////////////////////////////////////////////////////
	32
	33	ClassImp(AliITSVertexerZ)
	34
	35
	36
	37	//______________________________________________________________________
8221b41b	38	AliITSVertexerZ::AliITSVertexerZ():AliITSVertexer(),
	39	fFirstL1(0),
	40	fLastL1(0),
	41	fFirstL2(0),
	42	fLastL2(0),
	43	fDiffPhiMax(0),
	44	fX0(0.),
	45	fY0(0.),
	46	fZFound(0),
	47	fZsig(0.),
	48	fZCombc(0),
	49	fZCombv(0),
	50	fZCombf(0),
	51	fLowLim(0.),
	52	fHighLim(0.),
	53	fStepCoarse(0),
	54	fStepFine(0),
	55	fTolerance(0.),
	56	fMaxIter(0){
0c6af5c9	57	// Default constructor
0c6af5c9	58	SetDiffPhiMax(0);
0c6af5c9	59	SetFirstLayerModules(0);
0c6af5c9	60	SetSecondLayerModules(0);
0c6af5c9	61	SetLowLimit(0.);
	62	SetHighLimit(0.);
	63	SetBinWidthCoarse(0.);
	64	SetBinWidthFine(0.);
	65	SetTolerance(0.);
ecc64c3f	66	SetPPsetting(0.,0.);
6fd990e3	67	ConfigIterations();
0c6af5c9	68	}
	69
	70	//______________________________________________________________________
8221b41b	71	AliITSVertexerZ::AliITSVertexerZ(TString fn, Float_t x0, Float_t y0):AliITSVertexer(fn),
	72	fFirstL1(0),
	73	fLastL1(0),
	74	fFirstL2(0),
	75	fLastL2(0),
	76	fDiffPhiMax(0),
	77	fX0(x0),
	78	fY0(y0),
	79	fZFound(0),
	80	fZsig(0.),
	81	fZCombc(0),
	82	fZCombv(0),
	83	fZCombf(0),
	84	fLowLim(0.),
	85	fHighLim(0.),
	86	fStepCoarse(0),
	87	fStepFine(0),
	88	fTolerance(0.),
	89	fMaxIter(0) {
0c6af5c9	90	// Standard Constructor
0c6af5c9	91	SetDiffPhiMax();
0c6af5c9	92	SetFirstLayerModules();
0c6af5c9	93	SetSecondLayerModules();
0c6af5c9	94	SetLowLimit();
	95	SetHighLimit();
	96	SetBinWidthCoarse();
	97	SetBinWidthFine();
	98	SetTolerance();
ecc64c3f	99	SetPPsetting();
6fd990e3	100	ConfigIterations();
0c6af5c9	101
	102	}
	103
	104	//______________________________________________________________________
8221b41b	105	AliITSVertexerZ::AliITSVertexerZ(const AliITSVertexerZ &vtxr) : AliITSVertexer(vtxr),
	106	fFirstL1(vtxr.fFirstL1),
	107	fLastL1(vtxr.fLastL1),
	108	fFirstL2(vtxr.fFirstL2),
	109	fLastL2(vtxr.fLastL2),
	110	fDiffPhiMax(vtxr.fDiffPhiMax),
	111	fX0(vtxr.fX0),
	112	fY0(vtxr.fY0),
	113	fZFound(vtxr.fZFound),
	114	fZsig(vtxr.fZsig),
	115	fZCombc(vtxr.fZCombc),
	116	fZCombv(vtxr.fZCombv),
	117	fZCombf(vtxr.fZCombf),
	118	fLowLim(vtxr.fLowLim),
	119	fHighLim(vtxr.fHighLim),
	120	fStepCoarse(vtxr.fStepCoarse),
	121	fStepFine(vtxr.fStepFine),
	122	fTolerance(vtxr.fTolerance),
	123	fMaxIter(vtxr.fMaxIter) {
0c6af5c9	124	// Copy constructor
8221b41b	125
0c6af5c9	126	}
	127
	128	//______________________________________________________________________
8221b41b	129	AliITSVertexerZ& AliITSVertexerZ::operator=(const AliITSVertexerZ& vtxr ){
0c6af5c9	130	// Assignment operator
8221b41b	131
	132	this->~AliITSVertexerZ();
	133	new(this) AliITSVertexerZ(vtxr);
0c6af5c9	134	return *this;
	135	}
	136
0c6af5c9	137	//______________________________________________________________________
0c6af5c9	138	AliITSVertexerZ::~AliITSVertexerZ() {
ecc64c3f	139	// Destructor
	140	delete fZCombc;
	141	delete fZCombf;
	142	delete fZCombv;
0c6af5c9	143	}
0c6af5c9	144
6fd990e3	145	//______________________________________________________________________
	146	void AliITSVertexerZ::ConfigIterations(Int_t noiter,Float_t *ptr){
	147	// configure the iterative procedure to gain efficiency for
	148	// pp events with very low multiplicity
	149	Float_t defaults[5]={0.05,0.1,0.2,0.3,0.5};
	150	fMaxIter=noiter;
	151	if(noiter>5){
	152	Error("ConfigIterations","Maximum number of iterations is 5\n");
	153	fMaxIter=5;
	154	}
	155	for(Int_t j=0;j<5;j++)fPhiDiffIter[j]=defaults[j];
	156	if(ptr)for(Int_t j=0;j<fMaxIter;j++)fPhiDiffIter[j]=ptr[j];
	157	}
	158
debb3dde	159	//______________________________________________________________________
	160	Int_t AliITSVertexerZ::GetPeakRegion(TH1F*h, Int_t &binmin, Int_t &binmax) const {
	161	// Finds a region around a peak in the Z histogram
	162	// Case of 2 peaks is treated
	163	Int_t imax=h->GetNbinsX();
	164	Float_t maxval=0;
	165	Int_t bi1=h->GetMaximumBin();
	166	Int_t bi2=0;
	167	for(Int_t i=imax;i>=1;i--){
	168	if(h->GetBinContent(i)>maxval){
	169	maxval=h->GetBinContent(i);
	170	bi2=i;
	171	}
	172	}
	173	Int_t npeaks=0;
	174
	175	if(bi1==bi2){
	176	binmin=bi1-3;
	177	binmax=bi1+3;
	178	npeaks=1;
	179	}else{
	180	TH1F copy = new TH1F(h);
	181	copy->SetBinContent(bi1,0.);
	182	copy->SetBinContent(bi2,0.);
	183	Int_t l1=TMath::Max(bi1-3,1);
	184	Int_t l2=TMath::Min(bi1+3,h->GetNbinsX());
	185	Float_t cont1=copy->Integral(l1,l2);
	186	Int_t ll1=TMath::Max(bi2-3,1);
	187	Int_t ll2=TMath::Min(bi2+3,h->GetNbinsX());
	188	Float_t cont2=copy->Integral(ll1,ll2);
	189	if(cont1>cont2){
	190	binmin=l1;
	191	binmax=l2;
	192	npeaks=1;
	193	}
	194	if(cont2>cont1){
	195	binmin=ll1;
	196	binmax=ll2;
	197	npeaks=1;
	198	}
	199	if(cont1==cont2){
	200	binmin=l1;
	201	binmax=ll2;
	202	if(bi2-bi1==1) npeaks=1;
	203	else npeaks=2;
	204	}
	205	delete copy;
	206	}
	207	return npeaks;
	208	}
0c6af5c9	209	//______________________________________________________________________
d681bb2d	210	AliESDVertex* AliITSVertexerZ::FindVertexForCurrentEvent(Int_t evnumber){
d681bb2d	211	// Defines the AliESDVertex for the current event
6fd990e3	212	VertexZFinder(evnumber);
	213	Int_t ntrackl=0;
	214	for(Int_t iteraz=0;iteraz<fMaxIter;iteraz++){
	215	if(fCurrentVertex) ntrackl=fCurrentVertex->GetNContributors();
	216	if(!fCurrentVertex \|\| ntrackl==0 \|\| ntrackl==-1){
	217	Float_t diffPhiMaxOrig=fDiffPhiMax;
	218	fDiffPhiMax=GetPhiMaxIter(iteraz);
	219	VertexZFinder(evnumber);
	220	fDiffPhiMax=diffPhiMaxOrig;
	221	}
	222	}
32e449be	223	FindMultiplicity(evnumber);
6fd990e3	224	return fCurrentVertex;
	225	}
	226
32e449be	227
	228
	229
6fd990e3	230	//______________________________________________________________________
	231	void AliITSVertexerZ::VertexZFinder(Int_t evnumber){
	232	// Defines the AliESDVertex for the current event
0c6af5c9	233	fCurrentVertex = 0;
0c6af5c9	234	AliRunLoader *rl =AliRunLoader::GetRunLoader();
7d62fb64	235	AliITSLoader* itsLoader = (AliITSLoader*)rl->GetLoader("ITSLoader");
32e449be	236	AliITSgeom* geom = itsLoader->GetITSgeom();
2257f27e	237	itsLoader->LoadRecPoints();
0c6af5c9	238	rl->GetEvent(evnumber);
0c6af5c9	239
7d62fb64	240	AliITSDetTypeRec detTypeRec;
0c6af5c9	241
0c6af5c9	242	TTree *tR = itsLoader->TreeR();
7d62fb64	243	detTypeRec.SetTreeAddressR(tR);
0c6af5c9	244	TClonesArray *itsRec = 0;
ecc64c3f	245	// lc and gc are local and global coordinates for layer 1
0c6af5c9	246	Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
0c6af5c9	247	Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
ecc64c3f	248	// lc2 and gc2 are local and global coordinates for layer 2
0c6af5c9	249	Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
	250	Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
	251
7d62fb64	252	itsRec = detTypeRec.RecPoints();
2257f27e	253	TBranch *branch;
00a7cc50	254	branch = tR->GetBranch("ITSRecPoints");
0c6af5c9	255
0c6af5c9	256	Int_t nrpL1 = 0;
0c6af5c9	257	Int_t nrpL2 = 0;
ecc64c3f	258	// By default fFirstL1=0 and fLastL1=79
ecc64c3f	259	// This loop counts the number of recpoints on layer1 (central modules)
0c6af5c9	260	for(Int_t module= fFirstL1; module<=fLastL1;module++){
ecc64c3f	261	// Keep only central modules
0c6af5c9	262	if(module%4==0 \|\| module%4==3)continue;
ecc64c3f	263	// cout<<"Procesing module "<<module<<" ";
0c6af5c9	264	branch->GetEvent(module);
ecc64c3f	265	// cout<<"Number of clusters "<<clusters->GetEntries()<<endl;
0c6af5c9	266	nrpL1+= itsRec->GetEntries();
7d62fb64	267	detTypeRec.ResetRecPoints();
0c6af5c9	268	}
ecc64c3f	269	//By default fFirstL2=80 and fLastL2=239
ecc64c3f	270	//This loop counts the number of RP on layer 2
0c6af5c9	271	for(Int_t module= fFirstL2; module<=fLastL2;module++){
	272	branch->GetEvent(module);
	273	nrpL2+= itsRec->GetEntries();
7d62fb64	274	detTypeRec.ResetRecPoints();
0c6af5c9	275	}
0c6af5c9	276	if(nrpL1 == 0 \|\| nrpL2 == 0){
0c6af5c9	277	ResetHistograms();
6fd990e3	278	itsLoader->UnloadRecPoints();
	279	fCurrentVertex = new AliESDVertex(0.,5.3,-2);
	280	return;
0c6af5c9	281	}
ecc64c3f	282	// The vertex finding is attempted only if the number of RP is !=0 on
	283	// both layers
	284	Float_t *xc1 = new Float_t [nrpL1]; // coordinates of the L1 Recpoints
0c6af5c9	285	Float_t *yc1 = new Float_t [nrpL1];
	286	Float_t *zc1 = new Float_t [nrpL1];
	287	Float_t *phi1 = new Float_t [nrpL1];
ecc64c3f	288	Float_t *xc2 = new Float_t [nrpL2]; // coordinates of the L1 Recpoints
0c6af5c9	289	Float_t *yc2 = new Float_t [nrpL2];
	290	Float_t *zc2 = new Float_t [nrpL2];
	291	Float_t *phi2 = new Float_t [nrpL2];
ecc64c3f	292	Int_t ind = 0;// running index for RP
	293	// Force a coarse bin size of 200 microns if the number of clusters on layer 2
	294	// is low
	295	if(nrpL2<fPPsetting[0])SetBinWidthCoarse(fPPsetting[1]);
	296	// By default nbinfine = (10+10)/0.0005=40000
	297	Int_t nbinfine = static_cast<Int_t>((fHighLim-fLowLim)/fStepFine);
	298	// By default nbincoarse=(10+10)/0.01=2000
	299	Int_t nbincoarse = static_cast<Int_t>((fHighLim-fLowLim)/fStepCoarse);
	300	// Set stepverycoarse = 3*fStepCoarse
	301	Int_t nbinvcoarse = static_cast<Int_t>((fHighLim-fLowLim)/(3.*fStepCoarse));
	302	if(fZCombc)delete fZCombc;
	303	fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);
	304	if(fZCombv)delete fZCombv;
	305	fZCombv = new TH1F("fZCombv","Z",nbinvcoarse,fLowLim,fLowLim+nbinvcoarse3.fStepCoarse);
	306	if(fZCombf)delete fZCombf;
	307	fZCombf = new TH1F("fZCombf","Z",nbinfine,fLowLim,fLowLim+nbinfine*fStepFine);
	308	// Loop on modules of layer 1 (restricted to central modules)
0c6af5c9	309	for(Int_t module= fFirstL1; module<=fLastL1;module++){
	310	if(module%4==0 \|\| module%4==3)continue;
	311	branch->GetEvent(module);
	312	Int_t nrecp1 = itsRec->GetEntries();
	313	for(Int_t j=0;j<nrecp1;j++){
	314	AliITSRecPoint recp = (AliITSRecPoint)itsRec->At(j);
ecc64c3f	315	// Local coordinates of this recpoint
00a7cc50	316	lc[0]=recp->GetDetLocalX();
00a7cc50	317	lc[2]=recp->GetDetLocalZ();
0c6af5c9	318	geom->LtoG(module,lc,gc);
ecc64c3f	319	// Global coordinates of this recpoints
	320	gc[0]-=fX0; // Possible beam offset in the bending plane
	321	gc[1]-=fY0; // " "
0c6af5c9	322	xc1[ind]=gc[0];
	323	yc1[ind]=gc[1];
	324	zc1[ind]=gc[2];
ecc64c3f	325	// azimuthal angle is computed in the interval 0 --> 2*pi
0c6af5c9	326	phi1[ind] = TMath::ATan2(gc[1],gc[0]);
	327	if(phi1[ind]<0)phi1[ind]=2*TMath::Pi()+phi1[ind];
	328	ind++;
	329	}
7d62fb64	330	detTypeRec.ResetRecPoints();
0c6af5c9	331	}
ecc64c3f	332	ind = 0; // the running index is reset for Layer 2
0c6af5c9	333	for(Int_t module= fFirstL2; module<=fLastL2;module++){
	334	branch->GetEvent(module);
	335	Int_t nrecp2 = itsRec->GetEntries();
	336	for(Int_t j=0;j<nrecp2;j++){
	337	AliITSRecPoint recp = (AliITSRecPoint)itsRec->At(j);
00a7cc50	338	lc[0]=recp->GetDetLocalX();
00a7cc50	339	lc[2]=recp->GetDetLocalZ();
0c6af5c9	340	geom->LtoG(module,lc,gc);
	341	gc[0]-=fX0;
	342	gc[1]-=fY0;
	343	xc2[ind]=gc[0];
	344	yc2[ind]=gc[1];
	345	zc2[ind]=gc[2];
	346	phi2[ind] = TMath::ATan2(gc[1],gc[0]);
	347	if(phi2[ind]<0)phi2[ind]=2*TMath::Pi()+phi2[ind];
	348	ind++;
	349	}
7d62fb64	350	detTypeRec.ResetRecPoints();
0c6af5c9	351	}
32e449be	352
c61c02f4	353	// Int_t nolines=0;
ecc64c3f	354	for(Int_t i=0;i<nrpL1;i++){ // loop on L1 RP
	355	Float_t r1=TMath::Sqrt(xc1[i]xc1[i]+yc1[i]yc1[i]); // radius L1 RP
	356	for(Int_t j=0;j<nrpL2;j++){ // loop on L2 RP
	357	Float_t diff = TMath::Abs(phi2[j]-phi1[i]); // diff in azimuth
	358	if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff; //diff<pi
	359	if(diff<fDiffPhiMax){ // cut on 10 milliradians by def.
	360	Float_t r2=TMath::Sqrt(xc2[j]xc2[j]+yc2[j]yc2[j]); // radius L2 RP
	361	Float_t zr0=(r2zc1[i]-r1zc2[j])/(r2-r1); //Z @ null radius
0c6af5c9	362	fZCombf->Fill(zr0);
0c6af5c9	363	fZCombc->Fill(zr0);
ecc64c3f	364	fZCombv->Fill(zr0);
32e449be	365	Double_t pA[3];
	366	Double_t pB[3];
	367	pA[0]=xc1[i];
	368	pA[1]=yc1[i];
	369	pA[2]=zc1[i];
	370	pB[0]=xc2[j];
	371	pB[1]=yc2[j];
	372	pB[2]=zc2[j];
c61c02f4	373	// MakeTracklet(pA,pB,nolines);
0c6af5c9	374	}
	375	}
	376	}
	377	delete [] xc1;
	378	delete [] yc1;
	379	delete [] zc1;
	380	delete [] phi1;
	381	delete [] xc2;
	382	delete [] yc2;
	383	delete [] zc2;
	384	delete [] phi2;
6fd990e3	385	Double_t contents = fZCombc->GetEntries()- fZCombc->GetBinContent(0)-fZCombc->GetBinContent(nbincoarse+1);
	386	if(contents<1.){
	387	// Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
0c6af5c9	388	ResetHistograms();
6fd990e3	389	itsLoader->UnloadRecPoints();
	390	fCurrentVertex = new AliESDVertex(0.,5.3,-1);
	391	return;
0c6af5c9	392	}
ecc64c3f	393
ecc64c3f	394	TH1F *hc = fZCombc;
ecc64c3f	395
ecc64c3f	396
debb3dde	397	if(hc->GetBinContent(hc->GetMaximumBin())<3)hc = fZCombv;
	398	Int_t binmin,binmax;
	399	Int_t nPeaks=GetPeakRegion(hc,binmin,binmax);
	400	if(nPeaks==2)AliWarning("2 peaks found");
	401	Int_t ii1=1+static_cast<Int_t>((hc->GetBinLowEdge(binmin)-fLowLim)/fStepFine);
	402	Int_t ii2=1+static_cast<Int_t>((hc->GetBinLowEdge(binmax)+hc->GetBinWidth(binmax)-fLowLim)/fStepFine);
	403	Float_t centre = 0.;
ecc64c3f	404	Int_t nn=0;
debb3dde	405	for(Int_t ii=ii1;ii<=ii2;ii++){
d531b7a1	406	centre+= fZCombf->GetBinCenter(ii)*fZCombf->GetBinContent(ii);
d531b7a1	407	nn+=static_cast<Int_t>(fZCombf->GetBinContent(ii));
ecc64c3f	408	}
3eb8b116	409	if (nn) centre/=nn;
debb3dde	410
ecc64c3f	411	// n1 is the bin number of fine histogram containing the point located 1 coarse bin less than "centre"
debb3dde	412	Int_t n1 = 1+static_cast<Int_t>((centre-hc->GetBinWidth(1)-fLowLim)/fStepFine);
ecc64c3f	413	// n2 is the bin number of fine histogram containing the point located 1 coarse bin more than "centre"
debb3dde	414	Int_t n2 = 1+static_cast<Int_t>((centre+hc->GetBinWidth(1)-fLowLim)/fStepFine);
ecc64c3f	415	if(n1<1)n1=1;
ecc64c3f	416	if(n2>nbinfine)n2=nbinfine;
0c6af5c9	417	Int_t niter = 0;
debb3dde	418	Bool_t goon = kTRUE;
ecc64c3f	419	Int_t num=0;
	420	Bool_t last = kFALSE;
	421
	422	while(goon \|\| last){
0c6af5c9	423	fZFound = 0.;
	424	fZsig = 0.;
	425	num=0;
ecc64c3f	426	// at the end of the loop:
	427	// fZFound = N*(Average Z) - where N is the number of tracklets
	428	// num=N
	429	// fZsig = NQ - where Q is the average of ZZ
0c6af5c9	430	for(Int_t n=n1;n<=n2;n++){
	431	fZFound+=fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
	432	num+=static_cast<Int_t>(fZCombf->GetBinContent(n));
	433	fZsig+=fZCombf->GetBinCenter(n)fZCombf->GetBinCenter(n)fZCombf->GetBinContent(n);
	434	}
	435	if(num<2){
12e26c26	436	fZsig = 5.3; // Default error from the beam sigmoid
0c6af5c9	437	}
ecc64c3f	438	else{
979e3647	439	Float_t radi = fZsig/(num-1)-fZFound*fZFound/num/(num-1);
ecc64c3f	440	// radi = square root of sample variance of Z
979e3647	441	if(radi>0.)fZsig=TMath::Sqrt(radi);
12e26c26	442	else fZsig=5.3; // Default error from the beam sigmoid
ecc64c3f	443	// fZfound - Average Z
0c6af5c9	444	fZFound/=num;
0c6af5c9	445	}
ecc64c3f	446	if(!last){
	447	// goon is true if the distance between the found Z and the lower bin differs from the distance between the found Z and
	448	// the upper bin by more than tolerance (0.002)
	449	goon = TMath::Abs(TMath::Abs(fZFound-fZCombf->GetBinCenter(n1))-TMath::Abs(fZFound-fZCombf->GetBinCenter(n2)))>fTolerance;
	450	// a window in the fine grained histogram is centered aroung the found Z. The width is 2 bins of
	451	// the coarse grained histogram
6fd990e3	452	if(num>0){
debb3dde	453	n1 = 1+static_cast<Int_t>((fZFound-hc->GetBinWidth(1)-fLowLim)/fStepFine);
6fd990e3	454	if(n1<1)n1=1;
debb3dde	455	n2 = 1+static_cast<Int_t>((fZFound+hc->GetBinWidth(1)-fLowLim)/fStepFine);
6fd990e3	456	if(n2>nbinfine)n2=nbinfine;
6fd990e3	457	}
6fd990e3	458	else {
debb3dde	459	n1 = 1+static_cast<Int_t>((centre-(niter+2)*hc->GetBinWidth(1)-fLowLim)/fStepFine);
debb3dde	460	n2 = 1+static_cast<Int_t>((centre+(niter+2)*hc->GetBinWidth(1)-fLowLim)/fStepFine);
6fd990e3	461	if(n1<1)n1=1;
6fd990e3	462	if(n2>nbinfine)n2=nbinfine;
ecc64c3f	463	}
ecc64c3f	464	niter++;
ecc64c3f	465	// no more than 10 adjusting iterations
debb3dde	466	if(niter>=10)goon = kFALSE;
ecc64c3f	467
	468	if((fZsig> 0.0001) && !goon && num>5){
	469	last = kTRUE;
debb3dde	470	n1 = 1+static_cast<Int_t>((fZFound-fZsig-fLowLim)/fStepFine);
ecc64c3f	471	if(n1<1)n1=1;
debb3dde	472	n2 = 1+static_cast<Int_t>((fZFound+fZsig-fLowLim)/fStepFine);
ecc64c3f	473	if(n2>nbinfine)n2=nbinfine;
ecc64c3f	474	}
	475	}
	476	else {
	477	last = kFALSE;
0c6af5c9	478	}
ecc64c3f	479
0c6af5c9	480	}
ecc64c3f	481	// if(fDebug>0)cout<<"Numer of Iterations "<<niter<<endl<<endl;
9679d5db	482	// if(num!=0)fZsig/=TMath::Sqrt(num);
12e26c26	483	if (fZsig<=0) fZsig=5.3; // Default error from the beam sigmoid
d681bb2d	484	fCurrentVertex = new AliESDVertex(fZFound,fZsig,num);
0c6af5c9	485	fCurrentVertex->SetTitle("vertexer: B");
0c6af5c9	486	ResetHistograms();
6fd990e3	487	itsLoader->UnloadRecPoints();
6fd990e3	488	return;
0c6af5c9	489	}
	490
	491	//_____________________________________________________________________
	492	void AliITSVertexerZ::ResetHistograms(){
	493	// delete TH1 data members
	494	if(fZCombc)delete fZCombc;
	495	if(fZCombf)delete fZCombf;
ecc64c3f	496	if(fZCombv)delete fZCombv;
0c6af5c9	497	fZCombc = 0;
0c6af5c9	498	fZCombf = 0;
ecc64c3f	499	fZCombv = 0;
0c6af5c9	500	}
	501
	502	//______________________________________________________________________
	503	void AliITSVertexerZ::FindVertices(){
	504	// computes the vertices of the events in the range FirstEvent - LastEvent
	505	AliRunLoader *rl = AliRunLoader::GetRunLoader();
	506	AliITSLoader* itsLoader = (AliITSLoader*) rl->GetLoader("ITSLoader");
	507	itsLoader->ReloadRecPoints();
	508	for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
6fd990e3	509	// cout<<"Processing event "<<i<<endl;
0c6af5c9	510	rl->GetEvent(i);
	511	FindVertexForCurrentEvent(i);
	512	if(fCurrentVertex){
	513	WriteCurrentVertex();
	514	}
ecc64c3f	515	/*
0c6af5c9	516	else {
	517	if(fDebug>0){
	518	cout<<"Vertex not found for event "<<i<<endl;
	519	cout<<"fZFound = "<<fZFound<<", fZsig= "<<fZsig<<endl;
	520	}
	521	}
ecc64c3f	522	*/
0c6af5c9	523	}
	524	}
	525
	526	//________________________________________________________
	527	void AliITSVertexerZ::PrintStatus() const {
	528	// Print current status
	529	cout <<"=======================================================\n";
	530	cout <<" First layer first and last modules: "<<fFirstL1<<", ";
	531	cout <<fLastL1<<endl;
	532	cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
	533	cout <<fLastL2<<endl;
	534	cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
	535	cout <<"Limits for Z histograms: "<<fLowLim<<"; "<<fHighLim<<endl;
	536	cout <<"Bin sizes for coarse and fine z histos "<<fStepCoarse<<"; "<<fStepFine<<endl;
	537	cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
	538	cout <<" Debug flag: "<<fDebug<<endl;
	539	cout <<"First event to be processed "<<fFirstEvent;
	540	cout <<"\n Last event to be processed "<<fLastEvent<<endl;
ecc64c3f	541	if(fZCombc){
	542	cout<<"fZCombc exists - entries="<<fZCombc->GetEntries()<<endl;
	543	}
	544	else{
	545	cout<<"fZCombc does not exist\n";
	546	}
	547	if(fZCombv){
	548	cout<<"fZCombv exists - entries="<<fZCombv->GetEntries()<<endl;
	549	}
	550	else{
	551	cout<<"fZCombv does not exist\n";
	552	}
	553	if(fZCombf){
	554	cout<<"fZCombf exists - entries="<<fZCombv->GetEntries()<<endl;
	555	}
	556	else{
	557	cout<<"fZCombf does not exist\n";
	558	}
0c6af5c9	559
	560	cout <<"=======================================================\n";
	561	}
	562