stdlib.h included (HP,Sun)
[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 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.42198*ZCentroid+7.44718e-01*TMath::Power(ZCentroid,2)
110                   -5.73426e-01*TMath::Power(ZCentroid,3)+2.01500e-01*TMath::Power(ZCentroid,4)          
111               -3.34118e-02*TMath::Power(ZCentroid,5)+2.20816e-03*TMath::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.1*TMath::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-02*Vzero[2]+1.49626e-02*TMath::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';        
228    cout << "\nCentroid and RMS of Zv distribution: \n";
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"); 
252
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