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

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

#include "AliRun.h"
#include "AliITS.h"
#include "AliITSgeom.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>


ClassImp(AliITSVertex)

//////////////////////////////////////////////////////////////////////
// AliITSVertex is a class for full 3D primary vertex finding       //
//                                                                  //                                                           //
// Version: 1                                                       //                                   //     
// Written by Giuseppe Lo Re and Francesco Riggi                    //
// Giuseppe.Lore@ct.infn.it                                         //                           //
// Franco.Riggi@ct.infn.it                                          //                           //
// Release date: May 2001                                           //                                                                                           //
//                                                                  //                                                   //
//                                                                  //                                                   //
//////////////////////////////////////////////////////////////////////


//______________________________________________________________________
AliITSVertex::AliITSVertex() {
    // Default Constructor

    fPosition   = 0;
    fResolution = 0;
    fSNR        = 0;
}
//______________________________________________________________________
AliITSVertex::~AliITSVertex() {
    // Default Constructor

    if(fPosition){ delete [] fPosition; fPosition   = 0;}
    if(fResolution){ delete []fResolution; fResolution = 0;}
    if(fSNR){ delete [] fSNR;fSNR        = 0;}
}
//______________________________________________________________________
void AliITSVertex::Exec(){

   fPosition = new Double_t[3];
   fResolution = new Double_t[3];
   fSNR = new Double_t[3];
  
   AliITS* aliits =(AliITS *)gAlice->GetDetector("ITS");
   AliITSgeom *g2 = ((AliITS*)aliits)->GetITSgeom(); 
   
   TClonesArray  *recpoints = aliits->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);
  
   for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++) 
   {
           aliits->ResetRecPoints(); 
       gAlice->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) {         
                        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++) 
   {
        aliits->ResetRecPoints(); 
        gAlice->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");
     
   fSNR[2] = fz->GetParameter(0)/fz->GetParameter(3);
   if(fSNR[2]<0.) { 
     cout << "\nNegative Signal to noise ratio for z!!!" << endl;
     cout << "The algorithm cannot find the z vertex position." << endl;
     exit(123456789);
   }
   else
   {
     fPosition[2] = fz->GetParameter(1);
     if(fPosition[2]<0) 
       {
         fPosition[2]=fPosition[2]-TMath::Abs(fPosition[2])*1.11/10000;
       }
     else
       {
         fPosition[2]=fPosition[2]+TMath::Abs(fPosition[2])*1.11/10000;
       }
   }
   fResolution[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))-fPosition[2])
                <4*fResolution[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"); 

   fSNR[0] = fx->GetParameter(0)/fx->GetParameter(3);
   if(fSNR[0]<0.) { 
     cout << "\nNegative Signal to noise ratio for x!!!" << endl;
     cout << "The algorithm cannot find the x vertex position." << endl;
     exit(123456789);  
   }
   else
   {
   fPosition[0]=fx->GetParameter(1);
   fResolution[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"); 

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

   niter++;
   
   Vzero[0] = fPosition[0];
   Vzero[1] = fPosition[1];
   Vzero[2] = fPosition[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;
   
}
//______________________________________________________________________
Double_t AliITSVertex::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;
}

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