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