[u/mrichter/AliRoot.git] / ITS / AliITSVertexerIons.cxx

#include <AliITSVertexerIons.h>
#include "stdlib.h"
#include <TMath.h>
#include <TRandom.h>
#include <TObjArray.h>
#include <TROOT.h>
#include <TFile.h>
#include <TTree.h>
#include <Riostream.h>

#include "AliRun.h"
#include "AliITS.h"
#include "AliITSgeom.h"
#include "AliITSLoader.h"
#include "AliITSRecPoint.h"
#include "AliGenerator.h"
#include "AliMagF.h"

#include <TH1.h>
#include <TF1.h>
#include <TCanvas.h>
#include <AliESDVertex.h>
#include <TObjArray.h>
#include <TObject.h>
#include <AliITSVertexerPPZ.h>
//////////////////////////////////////////////////////////////////////
// AliITSVertexerIons is a class for full 3D primary vertex         //
// finding optimized for Ion-Ion interactions                       //
//                                                                  // 
//                                                                  //
//                                                                  //
//                                                                  //
// Written by Giuseppe Lo Re and Francesco Riggi                    //
//
// Giuseppe.Lore@ct.infn.it                                         //
// Franco.Riggi@ct.infn.it                                          //
//                                                                  //
// Release date: May 2001                                           //
//                                                                  //
//                                                                  //       
//////////////////////////////////////////////////////////////////////

ClassImp(AliITSVertexerIons)


  //______________________________________________________________________
  AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer() {
  // Default Constructor

  fITS = 0;
  SetNpThreshold();
}

//______________________________________________________________________
AliITSVertexerIons::AliITSVertexerIons(TString fn):AliITSVertexer(fn) {
  // Standard constructor
  
  fITS = 0;
  SetNpThreshold();
}


//______________________________________________________________________
AliITSVertexerIons::~AliITSVertexerIons() {
  // Default Destructor
  fITS = 0;
}

//______________________________________________________________________
AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
  // Defines the AliESDVertex for the current event
  fCurrentVertex = 0;
  Double_t position[3];
  Double_t resolution[3];
  Double_t snr[3];
  AliRunLoader *rl = AliRunLoader::GetRunLoader();
  AliITSLoader* itsloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
  if(!fITS) { 
    fITS =(AliITS *)gAlice->GetDetector("ITS");
    if(!fITS) {
      Error("FindVertexForCurrentEvent","AliITS object was not found");
      return fCurrentVertex;
    }
  }
  fITS->SetTreeAddress();
  AliITSgeom *g2 = fITS->GetITSgeom(); 
  TClonesArray  *recpoints = fITS->RecPoints();
  TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
  TBranch *branch;
  if(fUseV2Clusters){
    branch = itsloader->TreeR()->GetBranch("Clusters");
    branch->SetAddress(&clusters);
  }
  else {
    branch = itsloader->TreeR()->GetBranch("ITSRecPoints");
  }

  AliITSRecPoint *pnt;

  Int_t nopoints1 = 0;                                         
  Int_t nopoints2 = 0;
  Double_t vzero[3];
  Double_t aspar[2];
     
  //------------ Rough Vertex evaluation ---------------------------------   

  Int_t i,npoints,ipoint,j,k,max,binmax;
  Double_t zCentroid;
  Float_t l[3], p[3];

  Double_t mxpiu = 0;
  Double_t mxmeno = 0;
  Double_t mypiu = 0;
  Double_t mymeno = 0;
   
  TH1F *hITSz1         = new TH1F("hITSz1","",100,-14.35,14.35);
  TTree *tr =  itsloader->TreeR(); 
  for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++) 
    {
      fITS->ResetRecPoints();
      tr->GetEvent(i); 
      if(fUseV2Clusters){
	Clusters2RecPoints(clusters,i,recpoints);
      }
      npoints = recpoints->GetEntries();
      for (ipoint=0;ipoint<npoints;ipoint++) {
	pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
	l[0]=pnt->GetX();
	l[1]=0;
	l[2]=pnt->GetZ();
	g2->LtoG(i, l, p);
	if(i<80 && TMath::Abs(p[2])<14.35) {         
	  if(p[0]>0) mxpiu++; 
	  if(p[0]<0) mxmeno++; 
	  if(p[1]>0) mypiu++; 
	  if(p[1]<0) mymeno++; 
	  hITSz1->Fill(p[2]);       
	}
	if(i>=80 && TMath::Abs(p[2])<14.35) nopoints2++;
      }       
    }  
  
  nopoints1 = (Int_t)(hITSz1->GetEntries()); 
  if (nopoints1 == 0) {
    delete hITSz1;
    return fCurrentVertex;
  }         
  aspar[0] = (mxpiu-mxmeno)/(mxpiu+mxmeno);
  aspar[1] = (mypiu-mymeno)/(mypiu+mymeno); 
   
  vzero[0] = 5.24441*aspar[0]; 
  vzero[1] = 5.24441*aspar[1]; 

  zCentroid=  TMath::Abs(hITSz1->GetMean()); 
  vzero[2] = -5.31040e-02+1.42198*zCentroid+7.44718e-01*TMath::Power(zCentroid,2)
    -5.73426e-01*TMath::Power(zCentroid,3)+2.01500e-01*TMath::Power(zCentroid,4)          
    -3.34118e-02*TMath::Power(zCentroid,5)+2.20816e-03*TMath::Power(zCentroid,6);
   
  if(hITSz1->GetMean()<0) vzero[2] = -vzero[2];
   
  //cout << "\nXvzero: " << vzero[0] << " cm" << "";
  //cout << "\nYvzero: " << vzero[1] << " cm" << "";
  //cout << "\nZvzero: " << vzero[2] << " cm" << "\n";

  delete hITSz1;

  Double_t dphi,r,deltaZ,a,b;
  Int_t np1=0;
  Int_t np2=0;
  Int_t niter=0;
   
  Double_t deltaPhiZ = 0.08;
   
  Float_t mag[3];
  Float_t origin[3];
  for(Int_t okm=0;okm<3;okm++) origin[okm]=0;
  gAlice->Field()->Field(origin,mag);

  deltaPhiZ = deltaPhiZ*TMath::Abs(mag[2])/2;
  Double_t deltaPhiXY = 1.0;   

  //cout << "\ndeltaPhiZ: " << deltaPhiZ << " deg" << "\n";   
  //cout << "deltaPhiXY: " << deltaPhiXY << " deg" << "\n";   
   
  Double_t *z1, *z2, *y1, *y2, *x1, *x2, *phi1, *phi2, *r1, *r2;
  z1=new Double_t[nopoints1];
  z2=new Double_t[nopoints2];
  y1=new Double_t[nopoints1];
  y2=new Double_t[nopoints2];
  x1=new Double_t[nopoints1];
  x2=new Double_t[nopoints2];
  phi1=new Double_t[nopoints1];
  phi2=new Double_t[nopoints2];
  r1=new Double_t[nopoints1];
  r2=new Double_t[nopoints2];
        
  deltaZ = 4.91617e-01+2.67567e-02*vzero[2]+1.49626e-02*TMath::Power(vzero[2],2); 
  Float_t multFactorZ = 28000./(Float_t)nopoints1;
  Int_t nbin=(Int_t)((deltaZ/0.005)/multFactorZ);
  Int_t nbinxy=250;
  Int_t *vectorBinZ,*vectorBinXY;
  vectorBinZ=new Int_t[nbin];
  vectorBinXY=new Int_t[nbinxy];
  Float_t f1= 0;
  Float_t f2= 0;
  Double_t sigma,averagebg;
     
  TH1D *hITSZv         = new TH1D("hITSZv","",nbin,vzero[2]-deltaZ,vzero[2]+deltaZ);
  TH1D *hITSXv         = new TH1D("hITSXv","",nbinxy,-3,3);
  TH1D *hITSYv         = new TH1D("hITSYv","",nbinxy,-3,3);

  //cout << "deltaZeta: " << deltaZ << " cm" << "\n";   
   
   
 start:   

  hITSZv->Add(hITSZv,-1);
  hITSXv->Add(hITSXv,-1);
  hITSYv->Add(hITSYv,-1);

  np1=np2=0;


  for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++) 
    {
      fITS->ResetRecPoints(); 
      itsloader->TreeR()->GetEvent(i);
     if(fUseV2Clusters){
	Clusters2RecPoints(clusters,i,recpoints);
      }
      npoints = recpoints->GetEntries();
      for (ipoint=0;ipoint<npoints;ipoint++) {
               
	pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
	l[0]=pnt->GetX();
	l[1]=0;
	l[2]=pnt->GetZ();
	g2->LtoG(i, l, p);
              
	if(i<80 && TMath::Abs(p[2])<14.35) {
	  p[0]=p[0]-vzero[0];
	  p[1]=p[1]-vzero[1];
	  r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
	  y1[np1]=p[1];
	  x1[np1]=p[0];
	  z1[np1]=p[2]; 
	  r1[np1]=r;
	  phi1[np1]=PhiFunc(p);
	  np1++;
	}

	if(i>=80 &&  TMath::Abs(p[2])<14.35) { 
	  p[0]=p[0]-vzero[0];
	  p[1]=p[1]-vzero[1];
	  r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
	  y2[np2]=p[1];
	  x2[np2]=p[0];
	  z2[np2]=p[2];
	  r2[np2]=r;
	  phi2[np2]=PhiFunc(p);
	  np2++;
	}
                        
      }
    }

  //------------------ Correlation between rec points ----------------------
    
  
  if(np1<fNpThreshold) {
    /*    cout << "Warning: AliITSVertexerIons finder is not reliable for low multiplicity events. May be you have to use AliITSVertexerPPZ.\n";
	  position[0]=vzero[0];
	  position[1]=vzero[1];
	  position[2]=vzero[2];
	  resolution[0]=-123;
	  resolution[1]=-123;
	  resolution[2]=-123;
	  snr[0]=-123;
	  snr[1]=-123;
	  snr[2]=-123;
	  Char_t name[30];
	  sprintf(name,"Vertex");
	  fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
	  return fCurrentVertex;*/

    Warning("FindVertexForCurrentEvent","AliITSVertexerIons finder is not reliable for low multiplicity events. Switching to AliITSVertexerPPZ with default parameters...\n");
    Warning("FindVertexForCurrentEvent","N rec points = %d - Threshold is %d",np1,fNpThreshold);
    AliITSVertexerPPZ *dovert = new AliITSVertexerPPZ("null");
    fCurrentVertex =dovert->FindVertexForCurrentEvent(rl->GetEventNumber());
    delete dovert;
    return fCurrentVertex;
  }


  for(j=0; j<(np2)-1; j++) {
    for(k=0; k<(np1)-1; k++) { 
      dphi=TMath::Abs(phi1[k]-phi2[j]);
      if(dphi>180) dphi = 360-dphi;
      if(dphi<deltaPhiZ && TMath::Abs((z1[k]-(z2[j]-z1[k])/((r2[j]/r1[k])-1))-vzero[2])
	 <deltaZ) hITSZv->Fill(z1[k]-(z2[j]-z1[k])/((r2[j]/r1[k])-1));        
    }
  }
      
  //cout << "\nNumber of used pairs: \n";

  a = vzero[2]-deltaZ; 
  b = vzero[2]+deltaZ; 
  max=(Int_t) hITSZv->GetMaximum();
  binmax=hITSZv->GetMaximumBin();
  sigma=0;
  for(i=0;i<nbin;i++) vectorBinZ[i]=(Int_t)hITSZv->GetBinContent(i);
  for(i=0;i<10;i++) f1=f1+vectorBinZ[i]/10;
  for(i=nbin-10;i<nbin;i++) f2=f2+vectorBinZ[i]/10;
  averagebg=(f1+f2)/2;
  for(i=0;i<nbin;i++) {
    if(vectorBinZ[i]-averagebg>(max-averagebg)*0.4 && 
       vectorBinZ[i]-averagebg<(max-averagebg)*0.7) {
      sigma=hITSZv->GetBinCenter(binmax)-hITSZv->GetBinCenter(i);
      sigma=TMath::Abs(sigma);
      if(sigma==0) sigma=0.05;
    }
  }
  
  
  TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",a,b);  
  
  fz->SetParameter(0,max);
  if(niter==0) {Double_t temp = hITSZv->GetBinCenter(binmax); vzero[2]=temp;}
  fz->SetParameter(1,vzero[2]);
  fz->SetParameter(2,sigma); 
  fz->SetParameter(3,averagebg);  
  fz->SetParLimits(2,0,999);
  fz->SetParLimits(3,0,999);
  
  hITSZv->Fit("fz","RQ0");
  
  snr[2] = fz->GetParameter(0)/fz->GetParameter(3);
  if(snr[2]<0.) { 
    Error("FindVertexForCurrentEvent","\nNegative Signal to noise ratio for z!!!\n");
    Error("FindVertexForCurrentEvent","The algorithm cannot find the z vertex position.\n");
    //    exit(123456789);
    return fCurrentVertex;
  }
  else {
    position[2] = fz->GetParameter(1);
    if(position[2]<0) 
      {
	position[2]=position[2]-TMath::Abs(position[2])*1.11/10000;
      }
    else
      {
	position[2]=position[2]+TMath::Abs(position[2])*1.11/10000;
      }
  }
  resolution[2] = fz->GetParError(1);
  
  for(j=0; j<(np2)-1; j++) {
    for(k=0; k<(np1)-1; k++) { 
      dphi=TMath::Abs(phi1[k]-phi2[j]);
      if(dphi>180) dphi = 360-dphi;
      if(dphi>deltaPhiXY) continue;
      if(TMath::Abs((z1[k]-(z2[j]-z1[k])/((r2[j]/r1[k])-1))-position[2])
	 <4*resolution[2]) {
	hITSXv->Fill(vzero[0]+(x2[j]-((x2[j]-x1[k])/(y2[j]-y1[k]))*y2[j]));
	hITSYv->Fill(vzero[1]+(y2[j]-((y2[j]-y1[k])/(x2[j]-x1[k]))*x2[j]));
      }
    }
  }
  
  TF1 *fx = new TF1("fx","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",vzero[0]-0.5,vzero[0]+0.5);  
  
  max=(Int_t) hITSXv->GetMaximum();
  binmax=hITSXv->GetMaximumBin();
  sigma=0;
  f1=f2=0;
  for(i=0;i<nbinxy;i++) vectorBinXY[i]=(Int_t)hITSXv->GetBinContent(i);
  for(i=0;i<10;i++) f1=f1+vectorBinXY[i]/10;
  for(i=nbinxy-10;i<nbinxy;i++) f2=f2+vectorBinXY[i]/10;
  averagebg=(f1+f2)/2;
  for(i=0;i<nbinxy;i++) {
    if(vectorBinXY[i]-averagebg>(max-averagebg)*0.4 && 
       vectorBinXY[i]-averagebg<(max-averagebg)*0.7) {
      sigma=hITSXv->GetBinCenter(binmax)-hITSXv->GetBinCenter(i);
      sigma=TMath::Abs(sigma);
      if(sigma==0) sigma=0.05;
    }
  }
  
  
  fx->SetParameter(0,max);
  fx->SetParameter(1,vzero[0]);
  fx->SetParameter(2,sigma); 
  fx->SetParameter(3,averagebg);
  
  hITSXv->Fit("fx","RQ0"); 
  
  snr[0] = fx->GetParameter(0)/fx->GetParameter(3);
  
  if(snr[0]<0.) { 
    Error("FindVertexForCurrentEvent","\nNegative Signal to noise ratio for x!\n");
    Error("FindVertexForCurrentEvent","The algorithm cannot find the x vertex position\n");
    //   exit(123456789);  
    return fCurrentVertex;
  }
  else {
    position[0]=fx->GetParameter(1);
    resolution[0]=fx->GetParError(1);
  }
  
  TF1 *fy = new TF1("fy","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",vzero[1]-0.5,vzero[1]+0.5);  
  
  max=(Int_t) hITSYv->GetMaximum();
  binmax=hITSYv->GetMaximumBin();
  sigma=0;
  f1=f2=0;
  for(i=0;i<nbinxy;i++) vectorBinXY[i]=(Int_t)hITSYv->GetBinContent(i);
  for(i=0;i<10;i++) f1=f1+vectorBinXY[i]/10;
  for(i=nbinxy-10;i<nbinxy;i++) f2=f2+vectorBinXY[i]/10;
  averagebg=(f1+f2)/2;
  for(i=0;i<nbinxy;i++) {
    if(vectorBinXY[i]-averagebg>(max-averagebg)*0.4 && 
       vectorBinXY[i]-averagebg<(max-averagebg)*0.7) {
      sigma=hITSYv->GetBinCenter(binmax)-hITSYv->GetBinCenter(i);
      sigma=TMath::Abs(sigma);
      if(sigma==0) sigma=0.05;
    }
  }
  
  fy->SetParameter(0,max);
  fy->SetParameter(1,vzero[1]);
  fy->SetParameter(2,sigma); 
  fy->SetParameter(3,averagebg);
  
  hITSYv->Fit("fy","RQ0"); 
  
  snr[1] = fy->GetParameter(0)/fy->GetParameter(3);
  if(snr[1]<0.) {   
    Error("FindVertexForCurrentEvent","\nNegative Signal to noise ratio for y!\n");
    Error("FindVertexForCurrentEvent","The algorithm cannot find the y vertex position.\n");
    //    exit(123456789); 
    return fCurrentVertex;
  }
  else {
    position[1]=fy->GetParameter(1);
    resolution[1]=fy->GetParError(1);
  }
  
    
  /*cout << "iter = " << niter << " -> x = " << position[0] << endl;
    cout << "iter = " << niter << " -> y = " << position[1] << endl;
    cout << "iter = " << niter << " -> z = " << position[2] << endl;
    cout << "---\n";*/

  niter++;
   
  vzero[0] = position[0];
  vzero[1] = position[1];
  vzero[2] = position[2];

  if(niter<3) goto start;  // number of iterations
   

  cout<<"Number of iterations: "<<niter<<endl;
  cout << "\nNo. of Points on the two pixel layers: "<<np1<<" "<<np2<<endl; 
 
  delete [] z1;
  delete [] z2;
  delete [] y1;
  delete [] y2;
  delete [] x1;
  delete [] x2;
  delete [] r1;
  delete [] r2;
  delete [] phi1;
  delete [] phi2;
  delete [] vectorBinZ;
  delete [] vectorBinXY;
   
  delete hITSZv;
  delete hITSXv;
  delete hITSYv;
  Char_t name[30];
  //  sprintf(name,"Vertex_%d",evnumber);
  if(fDebug>0)Info("FindVertexForCurrentEvent","Vertex found for event %d",evnumber);
  sprintf(name,"Vertex");
  fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
  return fCurrentVertex;
}

//______________________________________________________________________
Double_t AliITSVertexerIons::PhiFunc(Float_t p[]) {
  Double_t phi=0;

  if(p[1]>0 && p[0]>0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578);
  if(p[1]>0 && p[0]<0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+180;
  if(p[1]<0 && p[0]<0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+180;
  if(p[1]<0 && p[0]>0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+360;
  return phi;
}

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

//________________________________________________________
void AliITSVertexerIons::PrintStatus() const {
  // Print current status
  cout <<"=======================================================\n";
  cout <<" Debug flag: "<<fDebug<<endl;
  cout<<"First event to be processed "<<fFirstEvent;
  cout<<"\n Last event to be processed "<<fLastEvent<<endl;
  if(fCurrentVertex)fCurrentVertex->PrintStatus();
}
Commit	Line	Data
c5f0f3c1	1	#include <AliITSVertexerIons.h>
	2	#include "stdlib.h"
	3	#include <TMath.h>
	4	#include <TRandom.h>
	5	#include <TObjArray.h>
	6	#include <TROOT.h>
	7	#include <TFile.h>
	8	#include <TTree.h>
	9	#include <Riostream.h>
	10
	11	#include "AliRun.h"
	12	#include "AliITS.h"
	13	#include "AliITSgeom.h"
88cb7938	14	#include "AliITSLoader.h"
c5f0f3c1	15	#include "AliITSRecPoint.h"
	16	#include "AliGenerator.h"
	17	#include "AliMagF.h"
	18
	19	#include <TH1.h>
	20	#include <TF1.h>
	21	#include <TCanvas.h>
d681bb2d	22	#include <AliESDVertex.h>
c5f0f3c1	23	#include <TObjArray.h>
c5f0f3c1	24	#include <TObject.h>
83028144	25	#include <AliITSVertexerPPZ.h>
c5f0f3c1	26	//////////////////////////////////////////////////////////////////////
	27	// AliITSVertexerIons is a class for full 3D primary vertex //
	28	// finding optimized for Ion-Ion interactions //
	29	// //
	30	// //
	31	// //
	32	// //
	33	// Written by Giuseppe Lo Re and Francesco Riggi //
83028144	34	//
c5f0f3c1	35	// Giuseppe.Lore@ct.infn.it //
c5f0f3c1	36	// Franco.Riggi@ct.infn.it //
	37	// //
	38	// Release date: May 2001 //
	39	// //
	40	// //
	41	//////////////////////////////////////////////////////////////////////
	42
	43	ClassImp(AliITSVertexerIons)
	44
	45
	46
83028144	47	//______________________________________________________________________
88cb7938	48	AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer() {
c5f0f3c1	49	// Default Constructor
	50
	51	fITS = 0;
83028144	52	SetNpThreshold();
c5f0f3c1	53	}
	54
	55	//______________________________________________________________________
88cb7938	56	AliITSVertexerIons::AliITSVertexerIons(TString fn):AliITSVertexer(fn) {
c5f0f3c1	57	// Standard constructor
88cb7938	58
c5f0f3c1	59	fITS = 0;
83028144	60	SetNpThreshold();
c5f0f3c1	61	}
	62
	63
	64	//______________________________________________________________________
	65	AliITSVertexerIons::~AliITSVertexerIons() {
	66	// Default Destructor
	67	fITS = 0;
	68	}
	69
	70	//______________________________________________________________________
d681bb2d	71	AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
d681bb2d	72	// Defines the AliESDVertex for the current event
c5f0f3c1	73	fCurrentVertex = 0;
83028144	74	Double_t position[3];
	75	Double_t resolution[3];
	76	Double_t snr[3];
88cb7938	77	AliRunLoader *rl = AliRunLoader::GetRunLoader();
83028144	78	AliITSLoader* itsloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
c5f0f3c1	79	if(!fITS) {
	80	fITS =(AliITS *)gAlice->GetDetector("ITS");
	81	if(!fITS) {
	82	Error("FindVertexForCurrentEvent","AliITS object was not found");
	83	return fCurrentVertex;
	84	}
	85	}
88cb7938	86	fITS->SetTreeAddress();
	87	AliITSgeom *g2 = fITS->GetITSgeom();
	88	TClonesArray *recpoints = fITS->RecPoints();
2257f27e	89	TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
	90	TBranch *branch;
	91	if(fUseV2Clusters){
	92	branch = itsloader->TreeR()->GetBranch("Clusters");
	93	branch->SetAddress(&clusters);
	94	}
	95	else {
	96	branch = itsloader->TreeR()->GetBranch("ITSRecPoints");
	97	}
	98
88cb7938	99	AliITSRecPoint *pnt;
88cb7938	100
83028144	101	Int_t nopoints1 = 0;
	102	Int_t nopoints2 = 0;
	103	Double_t vzero[3];
	104	Double_t aspar[2];
c5f0f3c1	105
88cb7938	106	//------------ Rough Vertex evaluation ---------------------------------
c5f0f3c1	107
83028144	108	Int_t i,npoints,ipoint,j,k,max,binmax;
83028144	109	Double_t zCentroid;
88cb7938	110	Float_t l[3], p[3];
c5f0f3c1	111
88cb7938	112	Double_t mxpiu = 0;
	113	Double_t mxmeno = 0;
	114	Double_t mypiu = 0;
	115	Double_t mymeno = 0;
c5f0f3c1	116
88cb7938	117	TH1F *hITSz1 = new TH1F("hITSz1","",100,-14.35,14.35);
83028144	118	TTree *tr = itsloader->TreeR();
88cb7938	119	for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++)
	120	{
	121	fITS->ResetRecPoints();
2257f27e	122	tr->GetEvent(i);
	123	if(fUseV2Clusters){
	124	Clusters2RecPoints(clusters,i,recpoints);
	125	}
88cb7938	126	npoints = recpoints->GetEntries();
	127	for (ipoint=0;ipoint<npoints;ipoint++) {
	128	pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
	129	l[0]=pnt->GetX();
	130	l[1]=0;
	131	l[2]=pnt->GetZ();
	132	g2->LtoG(i, l, p);
	133	if(i<80 && TMath::Abs(p[2])<14.35) {
	134	if(p[0]>0) mxpiu++;
	135	if(p[0]<0) mxmeno++;
	136	if(p[1]>0) mypiu++;
	137	if(p[1]<0) mymeno++;
	138	hITSz1->Fill(p[2]);
	139	}
83028144	140	if(i>=80 && TMath::Abs(p[2])<14.35) nopoints2++;
88cb7938	141	}
88cb7938	142	}
c5f0f3c1	143
83028144	144	nopoints1 = (Int_t)(hITSz1->GetEntries());
2257f27e	145	if (nopoints1 == 0) {
	146	delete hITSz1;
	147	return fCurrentVertex;
	148	}
83028144	149	aspar[0] = (mxpiu-mxmeno)/(mxpiu+mxmeno);
83028144	150	aspar[1] = (mypiu-mymeno)/(mypiu+mymeno);
c5f0f3c1	151
83028144	152	vzero[0] = 5.24441*aspar[0];
83028144	153	vzero[1] = 5.24441*aspar[1];
c5f0f3c1	154
83028144	155	zCentroid= TMath::Abs(hITSz1->GetMean());
	156	vzero[2] = -5.31040e-02+1.42198zCentroid+7.44718e-01TMath::Power(zCentroid,2)
	157	-5.73426e-01TMath::Power(zCentroid,3)+2.01500e-01TMath::Power(zCentroid,4)
	158	-3.34118e-02TMath::Power(zCentroid,5)+2.20816e-03TMath::Power(zCentroid,6);
c5f0f3c1	159
83028144	160	if(hITSz1->GetMean()<0) vzero[2] = -vzero[2];
c5f0f3c1	161
e83fe598	162	//cout << "\nXvzero: " << vzero[0] << " cm" << "";
	163	//cout << "\nYvzero: " << vzero[1] << " cm" << "";
	164	//cout << "\nZvzero: " << vzero[2] << " cm" << "\n";
c5f0f3c1	165
88cb7938	166	delete hITSz1;
c5f0f3c1	167
83028144	168	Double_t dphi,r,deltaZ,a,b;
88cb7938	169	Int_t np1=0;
	170	Int_t np2=0;
	171	Int_t niter=0;
c5f0f3c1	172
83028144	173	Double_t deltaPhiZ = 0.08;
c5f0f3c1	174
83028144	175	Float_t mag[3];
88cb7938	176	Float_t origin[3];
88cb7938	177	for(Int_t okm=0;okm<3;okm++) origin[okm]=0;
83028144	178	gAlice->Field()->Field(origin,mag);
c5f0f3c1	179
e83fe598	180	deltaPhiZ = deltaPhiZ*TMath::Abs(mag[2])/2;
83028144	181	Double_t deltaPhiXY = 1.0;
c5f0f3c1	182
e83fe598	183	//cout << "\ndeltaPhiZ: " << deltaPhiZ << " deg" << "\n";
e83fe598	184	//cout << "deltaPhiXY: " << deltaPhiXY << " deg" << "\n";
c5f0f3c1	185
83028144	186	Double_t z1, z2, y1, y2, x1, x2, phi1, phi2, r1, r2;
	187	z1=new Double_t[nopoints1];
	188	z2=new Double_t[nopoints2];
	189	y1=new Double_t[nopoints1];
	190	y2=new Double_t[nopoints2];
	191	x1=new Double_t[nopoints1];
	192	x2=new Double_t[nopoints2];
	193	phi1=new Double_t[nopoints1];
	194	phi2=new Double_t[nopoints2];
	195	r1=new Double_t[nopoints1];
	196	r2=new Double_t[nopoints2];
c5f0f3c1	197
83028144	198	deltaZ = 4.91617e-01+2.67567e-02vzero[2]+1.49626e-02TMath::Power(vzero[2],2);
	199	Float_t multFactorZ = 28000./(Float_t)nopoints1;
	200	Int_t nbin=(Int_t)((deltaZ/0.005)/multFactorZ);
88cb7938	201	Int_t nbinxy=250;
83028144	202	Int_t vectorBinZ,vectorBinXY;
	203	vectorBinZ=new Int_t[nbin];
	204	vectorBinXY=new Int_t[nbinxy];
88cb7938	205	Float_t f1= 0;
88cb7938	206	Float_t f2= 0;
83028144	207	Double_t sigma,averagebg;
c5f0f3c1	208
83028144	209	TH1D *hITSZv = new TH1D("hITSZv","",nbin,vzero[2]-deltaZ,vzero[2]+deltaZ);
88cb7938	210	TH1D *hITSXv = new TH1D("hITSXv","",nbinxy,-3,3);
88cb7938	211	TH1D *hITSYv = new TH1D("hITSYv","",nbinxy,-3,3);
c5f0f3c1	212
e83fe598	213	//cout << "deltaZeta: " << deltaZ << " cm" << "\n";
c5f0f3c1	214
c5f0f3c1	215
88cb7938	216	start:
c5f0f3c1	217
83028144	218	hITSZv->Add(hITSZv,-1);
	219	hITSXv->Add(hITSXv,-1);
	220	hITSYv->Add(hITSYv,-1);
c5f0f3c1	221
88cb7938	222	np1=np2=0;
c5f0f3c1	223
c5f0f3c1	224
88cb7938	225	for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++)
	226	{
	227	fITS->ResetRecPoints();
83028144	228	itsloader->TreeR()->GetEvent(i);
2257f27e	229	if(fUseV2Clusters){
	230	Clusters2RecPoints(clusters,i,recpoints);
	231	}
88cb7938	232	npoints = recpoints->GetEntries();
88cb7938	233	for (ipoint=0;ipoint<npoints;ipoint++) {
c5f0f3c1	234
88cb7938	235	pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
	236	l[0]=pnt->GetX();
	237	l[1]=0;
	238	l[2]=pnt->GetZ();
	239	g2->LtoG(i, l, p);
c5f0f3c1	240
88cb7938	241	if(i<80 && TMath::Abs(p[2])<14.35) {
83028144	242	p[0]=p[0]-vzero[0];
83028144	243	p[1]=p[1]-vzero[1];
88cb7938	244	r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
83028144	245	y1[np1]=p[1];
	246	x1[np1]=p[0];
	247	z1[np1]=p[2];
88cb7938	248	r1[np1]=r;
	249	phi1[np1]=PhiFunc(p);
	250	np1++;
	251	}
	252
	253	if(i>=80 && TMath::Abs(p[2])<14.35) {
83028144	254	p[0]=p[0]-vzero[0];
83028144	255	p[1]=p[1]-vzero[1];
88cb7938	256	r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
83028144	257	y2[np2]=p[1];
	258	x2[np2]=p[0];
	259	z2[np2]=p[2];
88cb7938	260	r2[np2]=r;
	261	phi2[np2]=PhiFunc(p);
	262	np2++;
	263	}
c5f0f3c1	264
88cb7938	265	}
88cb7938	266	}
c5f0f3c1	267
88cb7938	268	//------------------ Correlation between rec points ----------------------
c5f0f3c1	269
83028144	270
	271	if(np1<fNpThreshold) {
	272	/* cout << "Warning: AliITSVertexerIons finder is not reliable for low multiplicity events. May be you have to use AliITSVertexerPPZ.\n";
	273	position[0]=vzero[0];
	274	position[1]=vzero[1];
	275	position[2]=vzero[2];
	276	resolution[0]=-123;
	277	resolution[1]=-123;
	278	resolution[2]=-123;
	279	snr[0]=-123;
	280	snr[1]=-123;
	281	snr[2]=-123;
	282	Char_t name[30];
	283	sprintf(name,"Vertex");
d681bb2d	284	fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
83028144	285	return fCurrentVertex;*/
	286
	287	Warning("FindVertexForCurrentEvent","AliITSVertexerIons finder is not reliable for low multiplicity events. Switching to AliITSVertexerPPZ with default parameters...\n");
	288	Warning("FindVertexForCurrentEvent","N rec points = %d - Threshold is %d",np1,fNpThreshold);
2257f27e	289	AliITSVertexerPPZ *dovert = new AliITSVertexerPPZ("null");
83028144	290	fCurrentVertex =dovert->FindVertexForCurrentEvent(rl->GetEventNumber());
	291	delete dovert;
	292	return fCurrentVertex;
	293	}
	294
88cb7938	295
	296	for(j=0; j<(np2)-1; j++) {
	297	for(k=0; k<(np1)-1; k++) {
	298	dphi=TMath::Abs(phi1[k]-phi2[j]);
	299	if(dphi>180) dphi = 360-dphi;
83028144	300	if(dphi<deltaPhiZ && TMath::Abs((z1[k]-(z2[j]-z1[k])/((r2[j]/r1[k])-1))-vzero[2])
83028144	301	<deltaZ) hITSZv->Fill(z1[k]-(z2[j]-z1[k])/((r2[j]/r1[k])-1));
88cb7938	302	}
88cb7938	303	}
c5f0f3c1	304
88cb7938	305	//cout << "\nNumber of used pairs: \n";
83028144	306
	307	a = vzero[2]-deltaZ;
	308	b = vzero[2]+deltaZ;
88cb7938	309	max=(Int_t) hITSZv->GetMaximum();
83028144	310	binmax=hITSZv->GetMaximumBin();
88cb7938	311	sigma=0;
83028144	312	for(i=0;i<nbin;i++) vectorBinZ[i]=(Int_t)hITSZv->GetBinContent(i);
	313	for(i=0;i<10;i++) f1=f1+vectorBinZ[i]/10;
	314	for(i=nbin-10;i<nbin;i++) f2=f2+vectorBinZ[i]/10;
	315	averagebg=(f1+f2)/2;
88cb7938	316	for(i=0;i<nbin;i++) {
83028144	317	if(vectorBinZ[i]-averagebg>(max-averagebg)*0.4 &&
	318	vectorBinZ[i]-averagebg<(max-averagebg)*0.7) {
	319	sigma=hITSZv->GetBinCenter(binmax)-hITSZv->GetBinCenter(i);
88cb7938	320	sigma=TMath::Abs(sigma);
	321	if(sigma==0) sigma=0.05;
	322	}
	323	}
83028144	324
83028144	325
88cb7938	326	TF1 fz = new TF1 ("fz","([0]exp(-0.5((x-[1])/[2])((x-[1])/[2])))+[3]",a,b);
83028144	327
88cb7938	328	fz->SetParameter(0,max);
83028144	329	if(niter==0) {Double_t temp = hITSZv->GetBinCenter(binmax); vzero[2]=temp;}
83028144	330	fz->SetParameter(1,vzero[2]);
88cb7938	331	fz->SetParameter(2,sigma);
83028144	332	fz->SetParameter(3,averagebg);
88cb7938	333	fz->SetParLimits(2,0,999);
88cb7938	334	fz->SetParLimits(3,0,999);
83028144	335
	336	hITSZv->Fit("fz","RQ0");
	337
	338	snr[2] = fz->GetParameter(0)/fz->GetParameter(3);
	339	if(snr[2]<0.) {
	340	Error("FindVertexForCurrentEvent","\nNegative Signal to noise ratio for z!!!\n");
	341	Error("FindVertexForCurrentEvent","The algorithm cannot find the z vertex position.\n");
	342	// exit(123456789);
	343	return fCurrentVertex;
88cb7938	344	}
83028144	345	else {
	346	position[2] = fz->GetParameter(1);
	347	if(position[2]<0)
	348	{
	349	position[2]=position[2]-TMath::Abs(position[2])*1.11/10000;
	350	}
	351	else
	352	{
	353	position[2]=position[2]+TMath::Abs(position[2])*1.11/10000;
	354	}
	355	}
	356	resolution[2] = fz->GetParError(1);
	357
88cb7938	358	for(j=0; j<(np2)-1; j++) {
	359	for(k=0; k<(np1)-1; k++) {
	360	dphi=TMath::Abs(phi1[k]-phi2[j]);
	361	if(dphi>180) dphi = 360-dphi;
83028144	362	if(dphi>deltaPhiXY) continue;
	363	if(TMath::Abs((z1[k]-(z2[j]-z1[k])/((r2[j]/r1[k])-1))-position[2])
	364	<4*resolution[2]) {
	365	hITSXv->Fill(vzero[0]+(x2[j]-((x2[j]-x1[k])/(y2[j]-y1[k]))*y2[j]));
	366	hITSYv->Fill(vzero[1]+(y2[j]-((y2[j]-y1[k])/(x2[j]-x1[k]))*x2[j]));
c5f0f3c1	367	}
88cb7938	368	}
88cb7938	369	}
83028144	370
	371	TF1 fx = new TF1("fx","([0]exp(-0.5((x-[1])/[2])((x-[1])/[2])))+[3]",vzero[0]-0.5,vzero[0]+0.5);
	372
88cb7938	373	max=(Int_t) hITSXv->GetMaximum();
83028144	374	binmax=hITSXv->GetMaximumBin();
88cb7938	375	sigma=0;
88cb7938	376	f1=f2=0;
83028144	377	for(i=0;i<nbinxy;i++) vectorBinXY[i]=(Int_t)hITSXv->GetBinContent(i);
	378	for(i=0;i<10;i++) f1=f1+vectorBinXY[i]/10;
	379	for(i=nbinxy-10;i<nbinxy;i++) f2=f2+vectorBinXY[i]/10;
	380	averagebg=(f1+f2)/2;
88cb7938	381	for(i=0;i<nbinxy;i++) {
83028144	382	if(vectorBinXY[i]-averagebg>(max-averagebg)*0.4 &&
	383	vectorBinXY[i]-averagebg<(max-averagebg)*0.7) {
	384	sigma=hITSXv->GetBinCenter(binmax)-hITSXv->GetBinCenter(i);
88cb7938	385	sigma=TMath::Abs(sigma);
	386	if(sigma==0) sigma=0.05;
	387	}
	388	}
83028144	389
83028144	390
88cb7938	391	fx->SetParameter(0,max);
83028144	392	fx->SetParameter(1,vzero[0]);
88cb7938	393	fx->SetParameter(2,sigma);
83028144	394	fx->SetParameter(3,averagebg);
	395
	396	hITSXv->Fit("fx","RQ0");
	397
	398	snr[0] = fx->GetParameter(0)/fx->GetParameter(3);
	399
	400	if(snr[0]<0.) {
	401	Error("FindVertexForCurrentEvent","\nNegative Signal to noise ratio for x!\n");
	402	Error("FindVertexForCurrentEvent","The algorithm cannot find the x vertex position\n");
	403	// exit(123456789);
	404	return fCurrentVertex;
88cb7938	405	}
83028144	406	else {
	407	position[0]=fx->GetParameter(1);
	408	resolution[0]=fx->GetParError(1);
	409	}
	410
	411	TF1 fy = new TF1("fy","([0]exp(-0.5((x-[1])/[2])((x-[1])/[2])))+[3]",vzero[1]-0.5,vzero[1]+0.5);
	412
88cb7938	413	max=(Int_t) hITSYv->GetMaximum();
83028144	414	binmax=hITSYv->GetMaximumBin();
88cb7938	415	sigma=0;
88cb7938	416	f1=f2=0;
83028144	417	for(i=0;i<nbinxy;i++) vectorBinXY[i]=(Int_t)hITSYv->GetBinContent(i);
	418	for(i=0;i<10;i++) f1=f1+vectorBinXY[i]/10;
	419	for(i=nbinxy-10;i<nbinxy;i++) f2=f2+vectorBinXY[i]/10;
	420	averagebg=(f1+f2)/2;
88cb7938	421	for(i=0;i<nbinxy;i++) {
83028144	422	if(vectorBinXY[i]-averagebg>(max-averagebg)*0.4 &&
	423	vectorBinXY[i]-averagebg<(max-averagebg)*0.7) {
	424	sigma=hITSYv->GetBinCenter(binmax)-hITSYv->GetBinCenter(i);
88cb7938	425	sigma=TMath::Abs(sigma);
	426	if(sigma==0) sigma=0.05;
	427	}
	428	}
83028144	429
88cb7938	430	fy->SetParameter(0,max);
83028144	431	fy->SetParameter(1,vzero[1]);
88cb7938	432	fy->SetParameter(2,sigma);
83028144	433	fy->SetParameter(3,averagebg);
	434
	435	hITSYv->Fit("fy","RQ0");
	436
	437	snr[1] = fy->GetParameter(0)/fy->GetParameter(3);
	438	if(snr[1]<0.) {
	439	Error("FindVertexForCurrentEvent","\nNegative Signal to noise ratio for y!\n");
	440	Error("FindVertexForCurrentEvent","The algorithm cannot find the y vertex position.\n");
	441	// exit(123456789);
	442	return fCurrentVertex;
88cb7938	443	}
83028144	444	else {
	445	position[1]=fy->GetParameter(1);
	446	resolution[1]=fy->GetParError(1);
	447	}
	448
	449
	450	/*cout << "iter = " << niter << " -> x = " << position[0] << endl;
	451	cout << "iter = " << niter << " -> y = " << position[1] << endl;
	452	cout << "iter = " << niter << " -> z = " << position[2] << endl;
	453	cout << "---\n";*/
88cb7938	454
88cb7938	455	niter++;
c5f0f3c1	456
83028144	457	vzero[0] = position[0];
	458	vzero[1] = position[1];
	459	vzero[2] = position[2];
c5f0f3c1	460
88cb7938	461	if(niter<3) goto start; // number of iterations
c5f0f3c1	462
c5f0f3c1	463
88cb7938	464	cout<<"Number of iterations: "<<niter<<endl;
83028144	465	cout << "\nNo. of Points on the two pixel layers: "<<np1<<" "<<np2<<endl;
c5f0f3c1	466
83028144	467	delete [] z1;
	468	delete [] z2;
	469	delete [] y1;
	470	delete [] y2;
	471	delete [] x1;
	472	delete [] x2;
88cb7938	473	delete [] r1;
	474	delete [] r2;
	475	delete [] phi1;
	476	delete [] phi2;
83028144	477	delete [] vectorBinZ;
83028144	478	delete [] vectorBinXY;
c5f0f3c1	479
88cb7938	480	delete hITSZv;
	481	delete hITSXv;
	482	delete hITSYv;
	483	Char_t name[30];
	484	// sprintf(name,"Vertex_%d",evnumber);
088e0b8d	485	if(fDebug>0)Info("FindVertexForCurrentEvent","Vertex found for event %d",evnumber);
88cb7938	486	sprintf(name,"Vertex");
d681bb2d	487	fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
88cb7938	488	return fCurrentVertex;
c5f0f3c1	489	}
	490
	491	//______________________________________________________________________
	492	Double_t AliITSVertexerIons::PhiFunc(Float_t p[]) {
88cb7938	493	Double_t phi=0;
c5f0f3c1	494
88cb7938	495	if(p[1]>0 && p[0]>0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578);
	496	if(p[1]>0 && p[0]<0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+180;
	497	if(p[1]<0 && p[0]<0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+180;
	498	if(p[1]<0 && p[0]>0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+360;
	499	return phi;
c5f0f3c1	500	}
	501
	502	//______________________________________________________________________
	503	void AliITSVertexerIons::FindVertices(){
	504	// computes the vertices of the events in the range FirstEvent - LastEvent
88cb7938	505	AliRunLoader *rl = AliRunLoader::GetRunLoader();
83028144	506	AliITSLoader* itsloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
83028144	507	itsloader->LoadRecPoints("read");
ab6a6f40	508	for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
88cb7938	509	rl->GetEvent(i);
c5f0f3c1	510	FindVertexForCurrentEvent(i);
	511	if(fCurrentVertex){
	512	WriteCurrentVertex();
	513	}
	514	else {
	515	if(fDebug>0){
	516	cout<<"Vertex not found for event "<<i<<endl;
	517	}
	518	}
	519	}
c5f0f3c1	520	}
	521
	522	//________________________________________________________
	523	void AliITSVertexerIons::PrintStatus() const {
	524	// Print current status
88cb7938	525	cout <<"=======================================================\n";
88cb7938	526	cout <<" Debug flag: "<<fDebug<<endl;
c5f0f3c1	527	cout<<"First event to be processed "<<fFirstEvent;
	528	cout<<"\n Last event to be processed "<<fLastEvent<<endl;
	529	if(fCurrentVertex)fCurrentVertex->PrintStatus();
	530	}