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