78c99e0d0a0bb49df808a898773c0ae29ac71c48
[u/mrichter/AliRoot.git] / ITS / AliITSVertex.cxx
1 #include "stdlib.h"
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 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 //////////////////////////////////////////////////////////////////////
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
54    Int_t NoPoints1 = 0;                                         
55    Int_t NoPoints2 = 0;
56    Double_t Vzero[3];
57    Double_t AsPar[2];
58      
59 //------------ Rough Vertex evaluation ---------------------------------   
60
61    Int_t i,npoints,ipoint,j,k,max,BinMax;
62    Double_t ZCentroid;
63    Float_t l[3], p[3];
64
65    Double_t mxpiu = 0;
66    Double_t mxmeno = 0;
67    Double_t mypiu = 0;
68    Double_t mymeno = 0;
69    
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            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++;
92        }       
93    }  
94   
95    NoPoints1 = (Int_t)(hITSz1->GetEntries()); 
96            
97    AsPar[0] = (mxpiu-mxmeno)/(mxpiu+mxmeno);
98    AsPar[1] = (mypiu-mymeno)/(mypiu+mymeno); 
99    
100    Vzero[0] = 5.24441*AsPar[0]; 
101    Vzero[1] = 5.24441*AsPar[1]; 
102
103    ZCentroid=  TMath::Abs(hITSz1->GetMean()); 
104    Vzero[2] = -5.31040e-02+1.42198*ZCentroid+7.44718e-01*TMath::Power(ZCentroid,2)
105                   -5.73426e-01*TMath::Power(ZCentroid,3)+2.01500e-01*TMath::Power(ZCentroid,4)          
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    
110    cout << "\nXvzero: " << Vzero[0] << " cm" << "";
111    cout << "\nYvzero: " << Vzero[1] << " cm" << "";
112    cout << "\nZvzero: " << Vzero[2] << " cm" << "\n";
113
114    delete hITSz1;
115
116    Double_t dphi,r,DeltaZ,a,b;
117    Int_t np1=0;
118    Int_t np2=0;
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);
127
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    
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];
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);
160
161 //   cout << "DeltaZeta: " << DeltaZ << " cm" << "\n";   
162    
163    
164    start:   
165
166    hITSZv->Add(hITSZv,-1.);
167    hITSXv->Add(hITSXv,-1.);
168    hITSYv->Add(hITSYv,-1.);
169
170    np1=np2=0;
171
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    }
213
214 //------------------ Correlation between rec points ----------------------
215     
216    cout << "\nNo. of Points on the two pixel layers: "<<np1<<" "<<np2<<endl; 
217
218    for(j=0; j<(np2)-1; j++) {
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       }
245    }
246    
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;
254    
255    TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",a,b);  
256    
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);
264    
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);
272    }
273    else
274    {
275      fPosition[2] = fz->GetParameter(1);
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    }
285    fResolution[2] = fz->GetParError(1);
286    
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    }
348
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"); 
383
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;
423    
424 }
425
426
427 //______________________________________________________________________________
428
429 AliITSVertex::~AliITSVertex() { 
430    delete [] fPosition;
431    delete [] fResolution;
432    delete [] fSNR;
433 }
434
435 //______________________________________________________________________________
436
437 Double_t AliITSVertex::PhiFunc(Float_t p[]) {
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;
444 }
445
446 //______________________________________________________________________________
447
448