[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 primary vertex finding.
//
// Version: 0
// Written by Giuseppe Lo Re and Francesco Riggi
// Giuseppe.Lore@ct.infn.it
// Franco.Riggi@ct.infn.it
// Marh 2001
//
//
///////////////////////////////////////////////////////////////////////


//______________________________________________________________________________

AliITSVertex::AliITSVertex() {   

   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;
   TRandom rnd;

   Int_t NoPoints1 = 0;						
   Int_t NoPoints2 = 0;
   Double_t Vzero[3];
   Double_t AsPar[2];
     
//------------Rough Vertex------------------------------------------------------   

   Int_t i,npoints,ipoint,j,k;
   Double_t ZCentroid;
   Float_t l[3], p[3];

   TH1F *hITSx1pos      = new TH1F("hITSx1pos","",100,0.,4.2);
   TH1F *hITSx1neg      = new TH1F("hITSx1neg","",100,-4.2,0.);
   TH1F *hITSy1pos      = new TH1F("hITSy1pos","",100,0.,4.2);
   TH1F *hITSy1neg      = new TH1F("hITSy1neg","",100,-4.2,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) hITSx1pos->Fill(p[0]); 
		        if(p[0]<0) hITSx1neg->Fill(p[0]); 
		        if(p[1]>0) hITSy1pos->Fill(p[1]); 
		        if(p[1]<0) hITSy1neg->Fill(p[1]); 
		   hITSz1->Fill(p[2]);       
	       }
		   if(i>=80 && TMath::Abs(p[2])<14.35) (NoPoints2)++;
       }       
   }  
  
   NoPoints1 = (Int_t)(hITSz1->GetEntries()); 
   	   
   Double_t mxpiu = hITSx1pos->GetEntries();
   Double_t mxmeno = hITSx1neg->GetEntries();
   Double_t mypiu = hITSy1pos->GetEntries();
   Double_t mymeno = hITSy1neg->GetEntries();

   AsPar[0] = (mxpiu-mxmeno)/(mxpiu+mxmeno);
   AsPar[1] = (mypiu-mymeno)/(mypiu+mymeno); 
   
   Vzero[0] = 5.24434*AsPar[0]; 
   Vzero[1] = 5.24434*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 << "\nCentroid and RMS of hIITSz1: \n";
   //cout << hITSz1->GetMean() << "       "  << hITSz1->GetRMS()  <<endl;
   //cout << "\nAsPar[0]: " << AsPar[0];
   //cout << "\nAsPar[1]: " << AsPar[1];
   //cout << "\nAsPar[2]: " << AsPar[2];
   cout << "\nXvzero: " << Vzero[0] << " cm" << "";
   cout << "\nYvzero: " << Vzero[1] << " cm" << "";
   cout << "\nZvzero: " << Vzero[2] << " cm" << "\n";

   delete hITSx1pos;
   delete hITSx1neg;
   delete hITSy1pos;
   delete hITSy1neg;
   delete hITSz1;

//-----------------------Get points---------------------------------------------

   Double_t dphi,DeltaPhi,DeltaZ,r;
   Int_t np1=0;
   Int_t np2=0;

   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];
	
   for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++) 
   {
	aliits->ResetRecPoints(); 
	gAlice->TreeR()->GetEvent(i);
	   npoints = recpoints->GetEntries();
	   for (ipoint=0;ipoint<npoints;ipoint++) {
	       
	        //if(rnd.Integer(10)>2) continue;
		    pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
	        l[0]=pnt->GetX();
	        l[1]=0;
	        l[2]=pnt->GetZ();
	        g2->LtoG(i, l, p);
		    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));
	      
			if(i<80 && TMath::Abs(p[2])<14.35) { 
			Z1[np1]=p[2]; 
			Y1[np1]=p[1];
			X1[np1]=p[0];
			r1[np1]=r;
			phi1[np1]=PhiFunc(p);
			np1++;
			}

			if(i>=80 &&  TMath::Abs(p[2])<14.35) { 
			Z2[np2]=p[2];
			Y2[np2]=p[1];
			X2[np2]=p[0];
			r2[np2]=r;
			phi2[np2]=PhiFunc(p);
	 		np2++;
	 		}
			
	  }
   }

//-------------------Correlations-----------------------------------------------
    
   //cout << "\nPoints number on the first and second layer: "<<np1<<" "<<np2<<endl;

   DeltaPhi = 0.085-0.1*TMath::Sqrt(Vzero[0]*Vzero[0]+Vzero[1]*Vzero[1]); 
   if(DeltaPhi<0) {
      cout << "\nThe algorith can't find the vertex. " << endl;
	  exit(123456789);
   }
   
   Float_t B[3];
   Float_t origin[3];
   for(Int_t okm=0;okm<3;okm++) origin[okm]=0;
   gAlice->Field()->Field(origin,B);

   DeltaPhi = DeltaPhi*B[2]/2;

   cout << "\nDeltaPhi: " << DeltaPhi << " deg" << "\n";   

   DeltaZ =
   0.2+2.91617e-01+2.67567e-02*Vzero[2]+1.49626e-02*TMath::Power(Vzero[2],2); 
     
   cout << "DeltaZeta: " << DeltaZ << " cm" << "\n";     

   TH1F *hITSZv         = new TH1F("hITSZv","",DeltaZ/0.005,Vzero[2]-DeltaZ,Vzero[2]+DeltaZ);

   for(j=0; j<(np2)-1; j++) {
	 for(k=0; k<(np1)-1; k++) { 
		if(TMath::Abs(Z1[k]-Z2[j])>13.) continue;
		dphi=TMath::Abs(phi1[k]-phi2[j]);
		if(dphi>180) dphi = 360-dphi;
		if(dphi<DeltaPhi) { 
		  if(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';	
   cout << "\nCentroid and RMS of Zv distribution: \n";
   cout <<  hITSZv->GetMean() << "       " << hITSZv->GetRMS() << "\n"<< "\n";

   delete [] Z1;
   delete [] Z2;
   delete [] Y1;
   delete [] Y2;
   delete [] X1;
   delete [] X2;
   delete [] r1;
   delete [] r2;
   delete [] phi1;
   delete [] phi2;

   Double_t a = Vzero[2]-DeltaZ; 
   Double_t b = Vzero[2]+DeltaZ; 
   
   TF1 *f4 = new TF1 ("f4","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",a,b);  
   f4->SetParameter(0,700.);
   f4->SetParameter(1,Vzero[2]);
   f4->SetParameter(2,0.05); // !!
   f4->SetParameter(3,50.);
   
   hITSZv->Fit("f4","RME0"); 

   
   delete hITSZv;
  
   fSNR[2] = f4->GetParameter(0)/f4->GetParameter(3);
   if(fSNR[2]<1.5) { 
     cout << "\nSignal to noise ratio too small!!!" << endl;
     cout << "The algorithm can't find the z vertex position." << endl;
	 exit(123456789);
   }
   else
   {
     fPosition[2] = (f4->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] = f4->GetParError(1);
   
   AliGenerator *gener = gAlice->Generator();
   Float_t Vx,Vy,Vz;
   gener->GetOrigin(Vx,Vy,Vz);
 
   fPosition[0]=(Double_t)Vx;
   fPosition[1]=(Double_t)Vy;

   fResolution[0] = 0;
   fResolution[1] = 0;

   fSNR[0] = 0;
   fSNR[1] = 0;
   
}


//______________________________________________________________________________

AliITSVertex::~AliITSVertex() { 
   delete [] fPosition;
   delete [] fResolution;
   delete [] fSNR;
}

//______________________________________________________________________________


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
	27	////////////////////////////////////////////////////////////////////////
	28	// AliITSVertex is a class for primary vertex finding.
	29	//
	30	// Version: 0
	31	// Written by Giuseppe Lo Re and Francesco Riggi
	32	// Giuseppe.Lore@ct.infn.it
	33	// Franco.Riggi@ct.infn.it
	34	// Marh 2001
	35	//
	36	//
	37	///////////////////////////////////////////////////////////////////////
	38
	39
	40	//______________________________________________________________________________
	41
	42	AliITSVertex::AliITSVertex() {
	43
	44	fPosition = new Double_t[3];
	45	fResolution = new Double_t[3];
	46	fSNR = new Double_t[3];
	47
	48	AliITS* aliits =(AliITS *)gAlice->GetDetector("ITS");
	49	AliITSgeom g2 = ((AliITS)aliits)->GetITSgeom();
	50
	51	TClonesArray *recpoints = aliits->RecPoints();
	52	AliITSRecPoint *pnt;
	53	TRandom rnd;
	54
	55	Int_t NoPoints1 = 0;
	56	Int_t NoPoints2 = 0;
	57	Double_t Vzero[3];
	58	Double_t AsPar[2];
	59
	60	//------------Rough Vertex------------------------------------------------------
	61
	62	Int_t i,npoints,ipoint,j,k;
	63	Double_t ZCentroid;
	64	Float_t l[3], p[3];
	65
66	TH1F *hITSx1pos = new TH1F("hITSx1pos","",100,0.,4.2);
67	TH1F *hITSx1neg = new TH1F("hITSx1neg","",100,-4.2,0.);
68	TH1F *hITSy1pos = new TH1F("hITSy1pos","",100,0.,4.2);
69	TH1F *hITSy1neg = new TH1F("hITSy1neg","",100,-4.2,0.);
70	TH1F *hITSz1 = new TH1F("hITSz1","",100,-14.35,14.35);
71
72	for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++)
73	{
74	aliits->ResetRecPoints();
75	gAlice->TreeR()->GetEvent(i);
76
77	npoints = recpoints->GetEntries();
78	for (ipoint=0;ipoint<npoints;ipoint++) {
79	pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
80	l[0]=pnt->GetX();
81	l[1]=0;
82	l[2]=pnt->GetZ();
83	g2->LtoG(i, l, p);
84	if(i<80 && TMath::Abs(p[2])<14.35) {
85	if(p[0]>0) hITSx1pos->Fill(p[0]);
86	if(p[0]<0) hITSx1neg->Fill(p[0]);
87	if(p[1]>0) hITSy1pos->Fill(p[1]);
88	if(p[1]<0) hITSy1neg->Fill(p[1]);
89	hITSz1->Fill(p[2]);
90	}
91	if(i>=80 && TMath::Abs(p[2])<14.35) (NoPoints2)++;
92	}
93	}
94
95	NoPoints1 = (Int_t)(hITSz1->GetEntries());
96
97	Double_t mxpiu = hITSx1pos->GetEntries();
98	Double_t mxmeno = hITSx1neg->GetEntries();
99	Double_t mypiu = hITSy1pos->GetEntries();
100	Double_t mymeno = hITSy1neg->GetEntries();
101
102	AsPar[0] = (mxpiu-mxmeno)/(mxpiu+mxmeno);
103	AsPar[1] = (mypiu-mymeno)/(mypiu+mymeno);
104
105	Vzero[0] = 5.24434*AsPar[0];
106	Vzero[1] = 5.24434*AsPar[1];
107
108	ZCentroid= TMath::Abs(hITSz1->GetMean());
109	Vzero[2] = -5.31040e-02+1.42198ZCentroid+7.44718e-01TMath::Power(ZCentroid,2)
110	-5.73426e-01TMath::Power(ZCentroid,3)+2.01500e-01TMath::Power(ZCentroid,4)
111	-3.34118e-02TMath::Power(ZCentroid,5)+2.20816e-03TMath::Power(ZCentroid,6);
112
113	if(hITSz1->GetMean()<0) Vzero[2] = -Vzero[2];
114
115	//cout << "\nCentroid and RMS of hIITSz1: \n";
116	//cout << hITSz1->GetMean() << " " << hITSz1->GetRMS() <<endl;
117	//cout << "\nAsPar[0]: " << AsPar[0];
118	//cout << "\nAsPar[1]: " << AsPar[1];
119	//cout << "\nAsPar[2]: " << AsPar[2];
120	cout << "\nXvzero: " << Vzero[0] << " cm" << "";
121	cout << "\nYvzero: " << Vzero[1] << " cm" << "";
122	cout << "\nZvzero: " << Vzero[2] << " cm" << "\n";
123
124	delete hITSx1pos;
125	delete hITSx1neg;
126	delete hITSy1pos;
127	delete hITSy1neg;
128	delete hITSz1;
129
130	//-----------------------Get points---------------------------------------------
131
132	Double_t dphi,DeltaPhi,DeltaZ,r;
133	Int_t np1=0;
134	Int_t np2=0;
135
136	Double_t Z1, Z2, Y1, Y2, X1, X2, phi1, phi2, r1, r2;
137	Z1=new Double_t[NoPoints1];
138	Z2=new Double_t[NoPoints2];
139	Y1=new Double_t[NoPoints1];
140	Y2=new Double_t[NoPoints2];
141	X1=new Double_t[NoPoints1];
142	X2=new Double_t[NoPoints2];
143	phi1=new Double_t[NoPoints1];
144	phi2=new Double_t[NoPoints2];
145	r1=new Double_t[NoPoints1];
146	r2=new Double_t[NoPoints2];
147
148	for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++)
149	{
150	aliits->ResetRecPoints();
151	gAlice->TreeR()->GetEvent(i);
152	npoints = recpoints->GetEntries();
153	for (ipoint=0;ipoint<npoints;ipoint++) {
154
155	//if(rnd.Integer(10)>2) continue;
156	pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
157	l[0]=pnt->GetX();
158	l[1]=0;
159	l[2]=pnt->GetZ();
160	g2->LtoG(i, l, p);
161	p[0]=p[0]-Vzero[0];
162	p[1]=p[1]-Vzero[1];
163	r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
164
165	if(i<80 && TMath::Abs(p[2])<14.35) {
166	Z1[np1]=p[2];
167	Y1[np1]=p[1];
168	X1[np1]=p[0];
169	r1[np1]=r;
170	phi1[np1]=PhiFunc(p);
171	np1++;
172	}
173
174	if(i>=80 && TMath::Abs(p[2])<14.35) {
175	Z2[np2]=p[2];
176	Y2[np2]=p[1];
177	X2[np2]=p[0];
178	r2[np2]=r;
179	phi2[np2]=PhiFunc(p);
180	np2++;
181	}
182
183	}
184	}
185
186	//-------------------Correlations-----------------------------------------------
187
188	//cout << "\nPoints number on the first and second layer: "<<np1<<" "<<np2<<endl;
189
190	DeltaPhi = 0.085-0.1TMath::Sqrt(Vzero[0]Vzero[0]+Vzero[1]*Vzero[1]);
191	if(DeltaPhi<0) {
192	cout << "\nThe algorith can't find the vertex. " << endl;
193	exit(123456789);
194	}
195
196	Float_t B[3];
197	Float_t origin[3];
198	for(Int_t okm=0;okm<3;okm++) origin[okm]=0;
199	gAlice->Field()->Field(origin,B);
200
201	DeltaPhi = DeltaPhi*B[2]/2;
202
203	cout << "\nDeltaPhi: " << DeltaPhi << " deg" << "\n";
204
205	DeltaZ =
206	0.2+2.91617e-01+2.67567e-02Vzero[2]+1.49626e-02TMath::Power(Vzero[2],2);
207
208	cout << "DeltaZeta: " << DeltaZ << " cm" << "\n";
209
210	TH1F *hITSZv = new TH1F("hITSZv","",DeltaZ/0.005,Vzero[2]-DeltaZ,Vzero[2]+DeltaZ);
211
212	for(j=0; j<(np2)-1; j++) {
213	for(k=0; k<(np1)-1; k++) {
214	if(TMath::Abs(Z1[k]-Z2[j])>13.) continue;
215	dphi=TMath::Abs(phi1[k]-phi2[j]);
216	if(dphi>180) dphi = 360-dphi;
217	if(dphi<DeltaPhi) {
218	if(TMath::Abs((Z1[k]-(Z2[j]-Z1[k])/((r2[j]/r1[k])-1))-Vzero[2])<DeltaZ)
219	hITSZv->Fill(Z1[k]-(Z2[j]-Z1[k])/((r2[j]/r1[k])-1));
220	}
221	}
222	}
223
224
225
226	cout << "\nNumber of used pairs: \n";
227	cout << hITSZv->GetEntries() << '\n' << '\n';
ba2e68c4	228	cout << "\nCentroid and RMS of Zv distribution: \n";
504de69b	229	cout << hITSZv->GetMean() << " " << hITSZv->GetRMS() << "\n"<< "\n";
	230
	231	delete [] Z1;
	232	delete [] Z2;
	233	delete [] Y1;
	234	delete [] Y2;
	235	delete [] X1;
	236	delete [] X2;
	237	delete [] r1;
	238	delete [] r2;
	239	delete [] phi1;
	240	delete [] phi2;
	241
	242	Double_t a = Vzero[2]-DeltaZ;
	243	Double_t b = Vzero[2]+DeltaZ;
	244
	245	TF1 f4 = new TF1 ("f4","([0]exp(-0.5((x-[1])/[2])((x-[1])/[2])))+[3]",a,b);
	246	f4->SetParameter(0,700.);
	247	f4->SetParameter(1,Vzero[2]);
	248	f4->SetParameter(2,0.05); // !!
	249	f4->SetParameter(3,50.);
	250
	251	hITSZv->Fit("f4","RME0");
ba2e68c4	252
504de69b	253
	254	delete hITSZv;
	255
	256	fSNR[2] = f4->GetParameter(0)/f4->GetParameter(3);
	257	if(fSNR[2]<1.5) {
	258	cout << "\nSignal to noise ratio too small!!!" << endl;
	259	cout << "The algorithm can't find the z vertex position." << endl;
	260	exit(123456789);
	261	}
	262	else
	263	{
	264	fPosition[2] = (f4->GetParameter(1));
	265	if(fPosition[2]<0)
	266	{
	267	fPosition[2]=fPosition[2]-TMath::Abs(fPosition[2])*1.11/10000;
	268	}
	269	else
	270	{
	271	fPosition[2]=fPosition[2]+TMath::Abs(fPosition[2])*1.11/10000;
	272	}
	273	}
	274	fResolution[2] = f4->GetParError(1);
	275
	276	AliGenerator *gener = gAlice->Generator();
	277	Float_t Vx,Vy,Vz;
	278	gener->GetOrigin(Vx,Vy,Vz);
	279
	280	fPosition[0]=(Double_t)Vx;
	281	fPosition[1]=(Double_t)Vy;
	282
	283	fResolution[0] = 0;
	284	fResolution[1] = 0;
	285
	286	fSNR[0] = 0;
	287	fSNR[1] = 0;
	288
	289	}
	290
	291
	292	//______________________________________________________________________________
	293
	294	AliITSVertex::~AliITSVertex() {
	295	delete [] fPosition;
	296	delete [] fResolution;
	297	delete [] fSNR;
	298	}
	299
	300	//______________________________________________________________________________
	301
	302
	303	Double_t AliITSVertex::PhiFunc(Float_t p[]) {
	304	Double_t phi=0;
	305	if(p[1]>0 && p[0]>0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578);
	306	if(p[1]>0 && p[0]<0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+180;
	307	if(p[1]<0 && p[0]<0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+180;
	308	if(p[1]<0 && p[0]>0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+360;
	309	return phi;
	310	}
	311	//
	312	//______________________________________________________________________________
	313