12 #include "AliITSgeom.h"
13 #include "AliITSRecPoint.h"
14 #include "AliGenerator.h"
20 #include <AliITSVertex.h>
21 #include <TObjArray.h>
25 ClassImp(AliITSVertex)
27 //////////////////////////////////////////////////////////////////////
28 // AliITSVertex is a class for full 3D primary vertex finding //
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 // //
37 //////////////////////////////////////////////////////////////////////
40 //______________________________________________________________________
41 AliITSVertex::AliITSVertex() {
42 // Default Constructor
48 //______________________________________________________________________
49 AliITSVertex::~AliITSVertex() {
50 // Default Constructor
52 if(fPosition){ delete [] fPosition; fPosition = 0;}
53 if(fResolution){ delete []fResolution; fResolution = 0;}
54 if(fSNR){ delete [] fSNR;fSNR = 0;}
56 //______________________________________________________________________
57 void AliITSVertex::Exec(){
59 fPosition = new Double_t[3];
60 fResolution = new Double_t[3];
61 fSNR = new Double_t[3];
63 AliITS* aliits =(AliITS *)gAlice->GetDetector("ITS");
64 AliITSgeom *g2 = ((AliITS*)aliits)->GetITSgeom();
66 TClonesArray *recpoints = aliits->RecPoints();
74 //------------ Rough Vertex evaluation ---------------------------------
76 Int_t i,npoints,ipoint,j,k,max,BinMax;
85 TH1F *hITSz1 = new TH1F("hITSz1","",100,-14.35,14.35);
87 for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++)
89 aliits->ResetRecPoints();
90 gAlice->TreeR()->GetEvent(i);
91 npoints = recpoints->GetEntries();
93 for (ipoint=0;ipoint<npoints;ipoint++) {
94 pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
99 if(i<80 && TMath::Abs(p[2])<14.35) {
106 if(i>=80 && TMath::Abs(p[2])<14.35) NoPoints2++;
110 NoPoints1 = (Int_t)(hITSz1->GetEntries());
112 AsPar[0] = (mxpiu-mxmeno)/(mxpiu+mxmeno);
113 AsPar[1] = (mypiu-mymeno)/(mypiu+mymeno);
115 Vzero[0] = 5.24441*AsPar[0];
116 Vzero[1] = 5.24441*AsPar[1];
118 ZCentroid= TMath::Abs(hITSz1->GetMean());
119 Vzero[2] = -5.31040e-02+1.42198*ZCentroid+7.44718e-01*TMath::Power(ZCentroid,2)
120 -5.73426e-01*TMath::Power(ZCentroid,3)+2.01500e-01*TMath::Power(ZCentroid,4)
121 -3.34118e-02*TMath::Power(ZCentroid,5)+2.20816e-03*TMath::Power(ZCentroid,6);
123 if(hITSz1->GetMean()<0) Vzero[2] = -Vzero[2];
125 /*cout << "\nXvzero: " << Vzero[0] << " cm" << "";
126 cout << "\nYvzero: " << Vzero[1] << " cm" << "";
127 cout << "\nZvzero: " << Vzero[2] << " cm" << "\n";*/
131 Double_t dphi,r,DeltaZ,a,b;
136 Double_t DeltaPhiZ = 0.08;
140 for(Int_t okm=0;okm<3;okm++) origin[okm]=0;
141 gAlice->Field()->Field(origin,B);
143 DeltaPhiZ = DeltaPhiZ*B[2]/2;
144 Double_t DeltaPhiXY = 1.0;
146 // cout << "\nDeltaPhiZ: " << DeltaPhiZ << " deg" << "\n";
147 // cout << "DeltaPhiXY: " << DeltaPhiXY << " deg" << "\n";
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];
161 DeltaZ = 4.91617e-01+2.67567e-02*Vzero[2]+1.49626e-02*TMath::Power(Vzero[2],2);
162 Float_t MulFactorZ = 28000./(Float_t)NoPoints1;
163 Int_t nbin=(Int_t)((DeltaZ/0.005)/MulFactorZ);
165 Int_t *VectorBinZ,*VectorBinXY;
166 VectorBinZ=new Int_t[nbin];
167 VectorBinXY=new Int_t[nbinxy];
170 Double_t sigma,MediaFondo;
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);
176 // cout << "DeltaZeta: " << DeltaZ << " cm" << "\n";
181 hITSZv->Add(hITSZv,-1.);
182 hITSXv->Add(hITSXv,-1.);
183 hITSYv->Add(hITSYv,-1.);
189 for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++)
191 aliits->ResetRecPoints();
192 gAlice->TreeR()->GetEvent(i);
193 npoints = recpoints->GetEntries();
194 for (ipoint=0;ipoint<npoints;ipoint++) {
196 pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
202 if(i<80 && TMath::Abs(p[2])<14.35) {
205 r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
210 phi1[np1]=PhiFunc(p);
214 if(i>=80 && TMath::Abs(p[2])<14.35) {
217 r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
222 phi2[np2]=PhiFunc(p);
229 //------------------ Correlation between rec points ----------------------
231 //cout << "\nNo. of Points on the two pixel layers: "<<np1<<" "<<np2<<endl;
233 for(j=0; j<(np2)-1; j++) {
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));
242 //cout << "\nNumber of used pairs: \n";
243 //cout << hITSZv->GetEntries() << '\n' << '\n';
246 max=(Int_t) hITSZv->GetMaximum();
247 BinMax=hITSZv->GetMaximumBin();
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;
262 /*cout << "f1 " <<f1 <<endl;
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;
268 cout << "Fondo " << MediaFondo<< endl;*/
270 TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",a,b);
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);
280 hITSZv->Fit("fz","RMEQ0");
282 fSNR[2] = fz->GetParameter(0)/fz->GetParameter(3);
284 cout << "\nNegative Signal to noise ratio for z!!!" << endl;
285 cout << "The algorithm cannot find the z vertex position." << endl;
290 fPosition[2] = fz->GetParameter(1);
293 fPosition[2]=fPosition[2]-TMath::Abs(fPosition[2])*1.11/10000;
297 fPosition[2]=fPosition[2]+TMath::Abs(fPosition[2])*1.11/10000;
300 fResolution[2] = fz->GetParError(1);
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])
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]));
318 ("fx","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",Vzero[0]-0.5,Vzero[0]+0.5);
320 max=(Int_t) hITSXv->GetMaximum();
321 BinMax=hITSXv->GetMaximumBin();
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;
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;*/
345 fx->SetParameter(0,max);
346 fx->SetParameter(1,Vzero[0]);
347 fx->SetParameter(2,sigma);
348 fx->SetParameter(3,MediaFondo);
350 hITSXv->Fit("fx","RMEQ0");
352 fSNR[0] = fx->GetParameter(0)/fx->GetParameter(3);
354 cout << "\nNegative Signal to noise ratio for x!!!" << endl;
355 cout << "The algorithm cannot find the x vertex position." << endl;
360 fPosition[0]=fx->GetParameter(1);
361 fResolution[0]=fx->GetParError(1);
365 ("fy","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",Vzero[1]-0.5,Vzero[1]+0.5);
367 max=(Int_t) hITSYv->GetMaximum();
368 BinMax=hITSYv->GetMaximumBin();
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;
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;*/
392 fy->SetParameter(0,max);
393 fy->SetParameter(1,Vzero[1]);
394 fy->SetParameter(2,sigma);
395 fy->SetParameter(3,MediaFondo);
397 hITSYv->Fit("fy","RMEQ0");
399 fSNR[1] = fy->GetParameter(0)/fy->GetParameter(3);
401 cout << "\nNegative Signal to noise ratio for y!!!" << endl;
402 cout << "The algorithm cannot find the y vertex position." << endl;
407 fPosition[1]=fy->GetParameter(1);
408 fResolution[1]=fy->GetParError(1);
413 Vzero[0] = fPosition[0];
414 Vzero[1] = fPosition[1];
415 Vzero[2] = fPosition[2];
417 if(niter<3) goto start; // number of iterations
420 cout<<"Number of iterations: "<<niter<<endl;
432 delete [] VectorBinZ;
433 delete [] VectorBinXY;
440 //______________________________________________________________________
441 Double_t AliITSVertex::PhiFunc(Float_t p[]) {
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;
451 //______________________________________________________________________