846c685daaf72e6ee858a04bcde872f923b66d06
[u/mrichter/AliRoot.git] / ITS / AliITSVertex.cxx
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.42198*ZCentroid+7.44718e-01*TMath::Power(ZCentroid,2)
109                   -5.73426e-01*TMath::Power(ZCentroid,3)+2.01500e-01*TMath::Power(ZCentroid,4)          
110               -3.34118e-02*TMath::Power(ZCentroid,5)+2.20816e-03*TMath::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.1*TMath::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-02*Vzero[2]+1.49626e-02*TMath::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