Revised algorithm extended to very peripheral heavy-ion collisions (from G. Lo Re)
[u/mrichter/AliRoot.git] / ITS / AliITSVertexerIons.cxx
1 #include <AliITSVertexerIons.h>
2 #include "stdlib.h"
3 #include <TMath.h>
4 #include <TRandom.h>
5 #include <TObjArray.h>
6 #include <TROOT.h>
7 #include <TFile.h>
8 #include <TTree.h>
9 #include <Riostream.h>
10
11 #include "AliRun.h"
12 #include "AliITS.h"
13 #include "AliITSgeom.h"
14 #include "AliITSLoader.h"
15 #include "AliITSRecPoint.h"
16 #include "AliGenerator.h"
17 #include "AliMagF.h"
18
19 #include <TH1.h>
20 #include <TF1.h>
21 #include <TCanvas.h>
22 #include <AliESDVertex.h>
23 #include <TObjArray.h>
24 #include <TObject.h>
25 #include <AliITSVertexerPPZ.h>
26 //////////////////////////////////////////////////////////////////////
27 // AliITSVertexerIons is a class for full 3D primary vertex         //
28 // finding optimized for Ion-Ion interactions                       //
29 //                                                                  // 
30 //                                                                  //
31 //                                                                  //
32 //                                                                  //
33 // Written by Giuseppe Lo Re and Francesco Riggi                    //
34 // Giuseppe.Lore@ct.infn.it                                         //
35 // Franco.Riggi@ct.infn.it                                          //
36 //                                                                  //
37 // Release date: Mar 2004                                           //
38 //                                                                  //
39 //                                                                  //       
40 //////////////////////////////////////////////////////////////////////
41
42 ClassImp(AliITSVertexerIons)
43
44
45
46 //______________________________________________________________________
47   AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer() {
48   // Default Constructor
49
50   fITS = 0;
51   SetNpThreshold();
52   SetMaxDeltaPhi();
53   SetMaxDeltaZ();
54 }
55
56 //______________________________________________________________________
57 AliITSVertexerIons::AliITSVertexerIons(TString fn):AliITSVertexer(fn) {
58   // Standard constructor
59   
60   fITS = 0;
61   SetNpThreshold();
62   SetMaxDeltaPhi();
63   SetMaxDeltaZ();}
64
65
66 //______________________________________________________________________
67 AliITSVertexerIons::~AliITSVertexerIons() {
68   // Default Destructor
69   fITS = 0;
70 }
71
72 //______________________________________________________________________
73 AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
74
75   fCurrentVertex = 0;
76
77   AliRunLoader *rl = AliRunLoader::GetRunLoader();
78   AliITSLoader* itsloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
79   if(!fITS) {
80     fITS =(AliITS *)gAlice->GetDetector("ITS");
81     if(!fITS) {
82       Error("FindVertexForCurrentEvent","AliITS object was not found");
83       return fCurrentVertex;
84     }
85   }
86   fITS->SetTreeAddress();
87   AliITSgeom *g2 = fITS->GetITSgeom();
88   TClonesArray  *recpoints = fITS->RecPoints();
89   AliITSRecPoint *pnt;
90   TTree *tr =  itsloader->TreeR();
91
92   Int_t npoints=0;
93   Int_t nopoints1=40000;
94   Int_t nopoints2=40000;
95   Float_t l[3], p[3];
96
97   Double_t *z1, *z2, *y1, *y2, *x1, *x2, *phi1, *phi2, *r1, *r2;
98   z1=new Double_t[nopoints1];
99   z2=new Double_t[nopoints2];
100   y1=new Double_t[nopoints1];
101   y2=new Double_t[nopoints2];
102   x1=new Double_t[nopoints1];
103   x2=new Double_t[nopoints2];
104   phi1=new Double_t[nopoints1];
105   phi2=new Double_t[nopoints2];
106   r1=new Double_t[nopoints1];
107   r2=new Double_t[nopoints2];
108
109   Double_t mxpiu = 0;
110   Double_t mxmeno = 0;
111   Double_t mypiu = 0;
112   Double_t mymeno = 0;
113   Double_t r=0;
114
115   Int_t np1=0, np2=0;
116   for(Int_t i=g2->GetStartSPD();i<=g2->GetLastSPD();i++) {
117     fITS->ResetRecPoints();
118     tr->GetEvent(i);
119     npoints = recpoints->GetEntries();
120     for (Int_t ipoint=0;ipoint<npoints;ipoint++) {
121       pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
122       l[0]=pnt->GetX();
123       l[1]=0;
124       l[2]=pnt->GetZ();
125       g2->LtoG(i, l, p);
126       r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
127
128       if(i<80 && TMath::Abs(p[2])<14.35)  {
129                 y1[np1]=p[1];
130                 x1[np1]=p[0];
131                 z1[np1]=p[2];
132                 if(p[0]>0) mxpiu++;
133                 if(p[0]<0) mxmeno++;
134                 if(p[1]>0) mypiu++;
135                 if(p[1]<0) mymeno++;
136                 np1++;
137       }
138       if(i>=80 && TMath::Abs(p[2])<14.35) {
139                 y2[np2]=p[1];
140                 x2[np2]=p[0];
141                 z2[np2]=p[2];
142                 np2++;
143       }
144     }
145   }
146
147   if(np1<fNpThreshold) {
148     Warning("FindVertexForCurrentEvent","AliITSVertexerIons finder is not reliable for low multiplicity events. Switching to AliITSVertexerPPZ with default parameters...\n");
149     Warning("FindVertexForCurrentEvent","N rec points = %d - Threshold is %d",np1,fNpThreshold);
150     AliITSVertexerPPZ *dovert = new AliITSVertexerPPZ("default");
151     fCurrentVertex =dovert->FindVertexForCurrentEvent(rl->GetEventNumber());
152     delete dovert;
153     return fCurrentVertex;
154   }
155
156   if(!np1 || !np2) {
157     Error("FindVertexForCurrentEvent","No points in the pixer layers");
158     return fCurrentVertex;
159   }  
160   if(fDebug>0)   cout << "N. points layer 1 and 2 : " << np1 << " " << np2 << endl;
161
162   Double_t asparx = (mxpiu-mxmeno)/(mxpiu+mxmeno);
163   Double_t aspary = (mypiu-mymeno)/(mypiu+mymeno);
164
165   Double_t x0 = 2.86*asparx;
166   Double_t y0 = 2.86*aspary;
167
168   if(fDebug>0)   cout << "Rough xy vertex = " << x0 << " " << y0 << endl;  
169
170   for(Int_t i=0;i<np1;i++) {
171     x1[i]-=x0;
172     y1[i]-=y0;
173     PhiFunc(x1[i],y1[i],phi1[i]);
174     r1[i]=TMath::Sqrt(x1[i]*x1[i]+y1[i]*y1[i]);
175   }
176   for(Int_t i=0;i<np2;i++) {
177     x2[i]-=x0;
178     y2[i]-=y0;
179     PhiFunc(x2[i],y2[i],phi2[i]);
180     r2[i]=TMath::Sqrt(x2[i]*x2[i]+y2[i]*y2[i]);
181   }
182    
183   Int_t nbinxy=400;
184   Int_t nbinz=400;
185   TH1F *hxv=new TH1F("hxv","",nbinxy,-5,5);
186   TH1F *hyv=new TH1F("hyv","",nbinxy,-5,5);
187   TH1F *hzv=new TH1F("hzv","",nbinz,-15,15);
188
189   Double_t dphi;
190   for(Int_t j=0; j<np1; j++) {
191     for(Int_t k=0; k<np2; k++) {
192       dphi=TMath::Abs(phi2[k]-phi1[j]);
193       if(dphi>180) dphi = 360-dphi;
194       if(dphi>fMaxDeltaPhi) continue; 
195       hzv->Fill(z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1));
196      }
197   }
198
199   // Fitting ...
200   Double_t max;
201   Int_t bin_max;
202   Double_t max_center;
203   
204   max = hzv->GetMaximum();
205   bin_max=hzv->GetMaximumBin();
206   max_center=hzv->GetBinCenter(bin_max);
207   Double_t dxy=1.5;
208   
209   TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",max_center-dxy,max_center+dxy);  
210   fz->SetParameter(0,max);
211   fz->SetParameter(1,max_center);
212   fz->SetParameter(2,0.1);
213   fz->SetLineColor(kRed);
214   hzv->Fit("fz","RQ0");
215
216   if(TMath::Abs(x0)>0.4 || TMath::Abs(y0)>0.4) {
217     Warning("FindVertexForCurrentEvent","AliITSVertexerIonsGeneral finder is not reliable for events with x and y vertex coordinates too far from the origin of the reference system. Their values will be set to 0.\n");
218     Double_t position[3]={0,0,fz->GetParameter(1)};
219     Double_t resolution[3]={0,0,0};
220     Double_t snr[3]={0,0,0};
221     Char_t name[30];
222     if(fDebug>0)Info("FindVertexForCurrentEvent","Vertex found for event %d",evnumber);
223     sprintf(name,"Vertex");
224     fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
225     return fCurrentVertex;
226   }
227   
228   for(Int_t j=0; j<np1; j++) {
229     for(Int_t k=0; k<np2; k++) {
230       if(TMath::Abs((z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1))-fz->GetParameter(1))>fMaxDeltaZ) continue;
231       if(y2[k]==y1[j]) continue;
232       hxv->Fill(x0+(x2[k]-((x2[k]-x1[j])/(y2[k]-y1[j]))*y2[k]));
233       if(x2[k]==x1[j]) continue;
234       hyv->Fill(y0+(y2[k]-((y2[k]-y1[j])/(x2[k]-x1[j]))*x2[k]));
235     }
236   }
237     
238   TF1 *fx = new TF1 ("fx","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]+[4]*x+[5]*x*x",x0-dxy,x0+dxy);  
239   fx->SetParameter(0,100);
240   Double_t dist=0.3;
241   Double_t x_approx=FindMaxAround(x0,hxv,dist);
242   if(fDebug>0) cout << "x_approx = " << x_approx << endl;
243   fx->SetParameter(1,x_approx);
244   Double_t dif_centroid=0.07;
245   fx->SetParLimits(1,x_approx-dif_centroid,x_approx+dif_centroid);
246   fx->SetParameter(2,0.1);
247   fx->SetLineColor(kRed);
248   hxv->Fit("fx","RQW0"); 
249   
250   TF1 *fy = new TF1 ("fy","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]+[4]*x+[5]*x*x",y0-dxy,y0+dxy);  
251   fy->SetParameter(0,100);
252   Double_t y_approx=FindMaxAround(y0,hyv,dist);
253   if(fDebug>0) cout << "y_approx = " << y_approx << endl;
254   fy->SetParameter(1,y_approx);
255   fy->SetParLimits(1,y_approx-dif_centroid,y_approx+dif_centroid);
256   fy->SetParameter(2,0.1);
257   fy->SetLineColor(kRed);
258   hyv->Fit("fy","RQW0");
259   
260   delete [] z1;
261   delete [] z2;
262   delete [] y1;
263   delete [] y2;
264   delete [] x1;
265   delete [] x2;
266   delete [] r1;
267   delete [] r2;
268   delete [] phi1;
269   delete [] phi2;
270   delete hxv;
271   delete hyv;
272   delete hzv;
273   
274   Double_t position[3]={fx->GetParameter(1),fy->GetParameter(1),fz->GetParameter(1)};
275   Double_t resolution[3]={fx->GetParameter(2),fy->GetParameter(2),fz->GetParameter(2)};
276   Double_t snr[3]={0,0,0};
277   
278   Char_t name[30];
279   if(fDebug>0)Info("FindVertexForCurrentEvent","Vertex found for event %d",evnumber);
280   sprintf(name,"Vertex");
281   fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
282   return fCurrentVertex;
283 }
284
285 //______________________________________________________________________
286 void AliITSVertexerIons::PhiFunc(Double_t &x,Double_t &y,Double_t &phi) {
287   if(y>0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578);
288   if(y>0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180;
289   if(y<0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180;
290   if(y<0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+360;;
291 }
292
293 //______________________________________________________________________
294 void AliITSVertexerIons::FindVertices(){
295   // computes the vertices of the events in the range FirstEvent - LastEvent
296   AliRunLoader *rl = AliRunLoader::GetRunLoader();
297   AliITSLoader* itsloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
298   itsloader->LoadRecPoints("read");
299   for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
300     rl->GetEvent(i);
301     FindVertexForCurrentEvent(i);
302     if(fCurrentVertex){
303       WriteCurrentVertex();
304     }
305     else {
306       if(fDebug>0){
307         cout<<"Vertex not found for event "<<i<<endl;
308       }
309     }
310   }
311 }
312
313 //________________________________________________________
314 void AliITSVertexerIons::PrintStatus() const {
315   // Print current status
316   cout <<"=======================================================\n";
317   cout <<" Debug flag: "<<fDebug<<endl;
318   cout<<"First event to be processed "<<fFirstEvent;
319   cout<<"\n Last event to be processed "<<fLastEvent<<endl;
320   if(fCurrentVertex)fCurrentVertex->PrintStatus();
321 }
322
323 //________________________________________________________
324 Double_t AliITSVertexerIons::FindMaxAround(Double_t point, TH1F *h, Double_t distance) {
325   Int_t max_content=0;
326   Int_t max_bin=0;
327   for(Int_t i=0;i<h->GetNbinsX();i++) {
328     Int_t content=(Int_t)h->GetBinContent(i);
329     Double_t center=(Double_t)h->GetBinCenter(i);
330     if(fabs(center-point)>distance) continue;
331     if(content>max_content) {max_content=content;max_bin=i;}
332   }
333   Double_t max=h->GetBinCenter(max_bin);
334   return max;
335 }
336
337