[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 "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(TFile *infile, TFile *outfile):AliITSVertexer(infile,outfile) {
  // Standard constructor
  if(!fInFile)Fatal("AliITSVertexerIons","No inputfile has been defined!");
  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];
 
  if(!fITS) { 
    fITS =(AliITS *)gAlice->GetDetector("ITS");
    if(!fITS) {
      Error("FindVertexForCurrentEvent","AliITS object was not found");
      return fCurrentVertex;
    }
  }
   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);
  
   for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++) 
   {
           fITS->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++) 
   {
        fITS->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");
     
   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);
   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
  if(!fOutFile){
    Error("FindVertices","The output file is not available - nothing done");
    return;
  }
  if(!gAlice)gAlice = (AliRun*)fInFile->Get("gAlice");
  if(!gAlice){
    Error("FindVertices","The AliRun object is not available - nothing done");
    return;
  }
  TDirectory *curdir = gDirectory;
  fInFile->cd();
  for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
    gAlice->GetEvent(i);
    FindVertexForCurrentEvent(i);
    if(fCurrentVertex){
      WriteCurrentVertex();
    }
    else {
      if(fDebug>0){
	cout<<"Vertex not found for event "<<i<<endl;
      }
    }
  }
  curdir->cd();
}

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