[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 <AliITSVertex.h>
#include <TObjArray.h>
#include <TObject.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;
}

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


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

//______________________________________________________________________
AliITSVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
  // Defines the AliITSVertex 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();
  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);    
      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()); 
           
  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 B[3];
  Float_t origin[3];
  for(Int_t okm=0;okm<3;okm++) origin[okm]=0;
  gAlice->Field()->Field(origin,B);

  DeltaPhiZ = DeltaPhiZ*B[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 MulFactorZ = 28000./(Float_t)NoPoints1;
  Int_t nbin=(Int_t)((DeltaZ/0.005)/MulFactorZ);
  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,MediaFondo;
     
  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);
      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 ----------------------
    
  //cout << "\nNo. of Points on the two pixel layers: "<<np1<<" "<<np2<<endl; 

  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";
  //cout << hITSZv->GetEntries() << '\n' << '\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;
  MediaFondo=(f1+f2)/2;
  for(i=0;i<nbin;i++) {
    if(VectorBinZ[i]-MediaFondo>(max-MediaFondo)*0.4 && 
       VectorBinZ[i]-MediaFondo<(max-MediaFondo)*0.7) {
      sigma=hITSZv->GetBinCenter(BinMax)-hITSZv->GetBinCenter(i);
      sigma=TMath::Abs(sigma);
      if(sigma==0) sigma=0.05;
    }
  }
   
  /*cout << "f1 " <<f1 <<endl;
    cout << "f2 " <<f2 <<endl;
    cout << "GetMaximumBin " <<hITSZv->GetMaximumBin() <<endl;
    cout << "nbin " <<nbin <<endl;
    cout << "max " << hITSZv->GetBinContent(BinMax)<<endl;
    cout << "sigma " <<sigma <<endl;
    cout << "Fondo " <<   MediaFondo<< endl;*/
   
  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,MediaFondo);  
  fz->SetParLimits(2,0,999);
  fz->SetParLimits(3,0,999);
   
  hITSZv->Fit("fz","RMEQ0");
     
  SNR[2] = fz->GetParameter(0)/fz->GetParameter(3);
  if(SNR[2]<0.) { 
    cout << "\nNegative Signal to noise ratio for z!!!" << endl;
    cout << "The algorithm cannot find the z vertex position." << endl;
    exit(123456789);
  }
  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;
  MediaFondo=(f1+f2)/2;
  for(i=0;i<nbinxy;i++) {
    if(VectorBinXY[i]-MediaFondo>(max-MediaFondo)*0.4 && 
       VectorBinXY[i]-MediaFondo<(max-MediaFondo)*0.7) {
      sigma=hITSXv->GetBinCenter(BinMax)-hITSXv->GetBinCenter(i);
      sigma=TMath::Abs(sigma);
      if(sigma==0) sigma=0.05;
    }
  }
   
  /*cout << "f1 " <<f1 <<endl;
    cout << "f2 " <<f2 <<endl;
    cout << "GetMaximumBin " <<hITSXv->GetMaximumBin() <<endl;
    cout << "max " << hITSXv->GetBinContent(BinMax)<<endl;
    cout << "sigma " <<sigma <<endl;
    cout << "Fondo " <<   MediaFondo<< endl;
    cout << "nbinxy " <<nbinxy <<endl;*/
   
  fx->SetParameter(0,max);
  fx->SetParameter(1,Vzero[0]);
  fx->SetParameter(2,sigma); 
  fx->SetParameter(3,MediaFondo);

  hITSXv->Fit("fx","RMEQ0"); 

  SNR[0] = fx->GetParameter(0)/fx->GetParameter(3);
  if(SNR[0]<0.) { 
    cout << "\nNegative Signal to noise ratio for x!!!" << endl;
    cout << "The algorithm cannot find the x vertex position." << endl;
    exit(123456789);  
  }
  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;
  MediaFondo=(f1+f2)/2;
  for(i=0;i<nbinxy;i++) {
    if(VectorBinXY[i]-MediaFondo>(max-MediaFondo)*0.4 && 
       VectorBinXY[i]-MediaFondo<(max-MediaFondo)*0.7) {
      sigma=hITSYv->GetBinCenter(BinMax)-hITSYv->GetBinCenter(i);
      sigma=TMath::Abs(sigma);
      if(sigma==0) sigma=0.05;
    }
  }
   
  /*cout << "f1 " <<f1 <<endl;
    cout << "f2 " <<f2 <<endl;
    cout << "GetMaximumBin " <<hITSYv->GetMaximumBin() <<endl;
    cout << "max " << hITSYv->GetBinContent(BinMax)<<endl;
    cout << "sigma " <<sigma <<endl;
    cout << "Fondo " <<   MediaFondo<< endl;
    cout << "nbinxy " <<nbinxy <<endl;*/
   
  fy->SetParameter(0,max);
  fy->SetParameter(1,Vzero[1]);
  fy->SetParameter(2,sigma); 
  fy->SetParameter(3,MediaFondo);
   
  hITSYv->Fit("fy","RMEQ0"); 

  SNR[1] = fy->GetParameter(0)/fy->GetParameter(3);
  if(SNR[1]<0.) { 
    cout << "\nNegative Signal to noise ratio for y!!!" << endl;
    cout << "The algorithm cannot find the y vertex position." << endl;
    exit(123456789); 
  }
  else
    {
      Position[1]=fy->GetParameter(1);
      Resolution[1]=fy->GetParError(1);
    }

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