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

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