For Pythia with tune don't switch off MI in ConfigHeavyFlavor
[u/mrichter/AliRoot.git] / ITS / AliITSVertexerIons.cxx
1 #include "AliRunLoader.h"
2 #include "AliITSDetTypeRec.h"
3 #include "AliITSVertexerIons.h"
4 #include "AliITSVertexerZ.h"
5 #include "AliESDVertex.h"
6 #include "AliITSgeomTGeo.h"
7 #include "AliITSLoader.h"
8 #include "AliITSRecPoint.h"
9 #include "AliLog.h"
10 #include <TTree.h>
11 #include <TH1.h>
12 #include <TF1.h>
13
14 //////////////////////////////////////////////////////////////////////
15 // AliITSVertexerIons is a class for full 3D primary vertex         //
16 // finding optimized for Ion-Ion interactions                       //
17 //                                                                  // 
18 //                                                                  //
19 //                                                                  //
20 //                                                                  //
21 // Written by Giuseppe Lo Re and Francesco Riggi                    //
22 // Giuseppe.Lore@ct.infn.it                                         //
23 // Franco.Riggi@ct.infn.it                                          //
24 //                                                                  //
25 // Release date: Mar 2004                                           //
26 //                                                                  //
27 //                                                                  //       
28 //////////////////////////////////////////////////////////////////////
29
30
31 ClassImp(AliITSVertexerIons)
32   
33 //______________________________________________________________________
34   AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer(),
35 fNpThreshold(0),
36 fMaxDeltaPhi(0),
37 fMaxDeltaZ(0){
38   // Default Constructor
39
40   //fITS = 0;
41   SetNpThreshold();
42   SetMaxDeltaPhi();
43   SetMaxDeltaZ();
44 }
45
46 //______________________________________________________________________
47 AliITSVertexerIons::~AliITSVertexerIons() {
48   // Default Destructor
49   //fITS = 0;
50 }
51
52 //______________________________________________________________________
53 AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(TTree *itsClusterTree){ 
54 // Defines the AliESDVertex for the current event
55
56   fCurrentVertex = 0;
57
58   AliITSDetTypeRec detTypeRec;
59
60   detTypeRec.SetTreeAddressR(itsClusterTree);
61
62   TClonesArray  *recpoints = detTypeRec.RecPoints();
63   AliITSRecPoint *pnt;
64
65
66   Int_t npoints=0;
67   Int_t nopoints1=40000;
68   Int_t nopoints2=40000;
69   Float_t p[3];
70
71   Double_t *z1, *z2, *y1, *y2, *x1, *x2, *phi1, *phi2, *r1, *r2;
72   z1=new Double_t[nopoints1];
73   z2=new Double_t[nopoints2];
74   y1=new Double_t[nopoints1];
75   y2=new Double_t[nopoints2];
76   x1=new Double_t[nopoints1];
77   x2=new Double_t[nopoints2];
78   phi1=new Double_t[nopoints1];
79   phi2=new Double_t[nopoints2];
80   r1=new Double_t[nopoints1];
81   r2=new Double_t[nopoints2];
82
83   Double_t mxpiu = 0;
84   Double_t mxmeno = 0;
85   Double_t mypiu = 0;
86   Double_t mymeno = 0;
87   Double_t r=0;
88
89   Int_t np1=0, np2=0;
90   for(Int_t i=AliITSgeomTGeo::GetModuleIndex(1,1,1);i<=AliITSgeomTGeo::GetModuleIndex(2,1,1)-1;i++) {
91     detTypeRec.ResetRecPoints();
92     itsClusterTree->GetEvent(i);
93     npoints = recpoints->GetEntries();
94     for (Int_t ipoint=0;ipoint<npoints;ipoint++) {
95       pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
96       /*
97       l[0]=pnt->GetDetLocalX();
98       l[1]=0;
99       l[2]=pnt->GetDetLocalZ();
100       g2->LtoG(i, l, p);
101       */
102       pnt->GetGlobalXYZ(p);
103       r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
104
105       if(i<80 && TMath::Abs(p[2])<14.35)  {
106                 y1[np1]=p[1];
107                 x1[np1]=p[0];
108                 z1[np1]=p[2];
109                 if(p[0]>0) mxpiu++;
110                 if(p[0]<0) mxmeno++;
111                 if(p[1]>0) mypiu++;
112                 if(p[1]<0) mymeno++;
113                 np1++;
114       }
115       if(i>=80 && TMath::Abs(p[2])<14.35) {
116                 y2[np2]=p[1];
117                 x2[np2]=p[0];
118                 z2[np2]=p[2];
119                 np2++;
120       }
121     }
122   }
123
124   if(np1<fNpThreshold) {
125     Warning("FindVertexForCurrentEvent","AliITSVertexerIons finder is not reliable for low multiplicity events. Switching to AliITSVertexerZ with default parameters...\n");
126     Warning("FindVertexForCurrentEvent","N rec points = %d - Threshold is %d",np1,fNpThreshold);
127     AliITSVertexerZ *dovert = new AliITSVertexerZ();
128     fCurrentVertex =dovert->FindVertexForCurrentEvent(itsClusterTree);
129     delete dovert;
130     return fCurrentVertex;
131   }
132
133   if(!np1 || !np2) {
134     Error("FindVertexForCurrentEvent","No points in the pixer layers");
135     return fCurrentVertex;
136   }  
137   AliDebug(1,Form("N. points layer 1 and 2 : %d %d",np1,np2));
138
139   Double_t asparx = (mxpiu-mxmeno)/(mxpiu+mxmeno);
140   Double_t aspary = (mypiu-mymeno)/(mypiu+mymeno);
141
142   Double_t x0 = 2.86*asparx;
143   Double_t y0 = 2.86*aspary;
144
145   AliDebug(1,Form("Rough xy vertex = %f %f",x0,y0));  
146
147   for(Int_t i=0;i<np1;i++) {
148     x1[i]-=x0;
149     y1[i]-=y0;
150     PhiFunc(x1[i],y1[i],phi1[i]);
151     r1[i]=TMath::Sqrt(x1[i]*x1[i]+y1[i]*y1[i]);
152   }
153   for(Int_t i=0;i<np2;i++) {
154     x2[i]-=x0;
155     y2[i]-=y0;
156     PhiFunc(x2[i],y2[i],phi2[i]);
157     r2[i]=TMath::Sqrt(x2[i]*x2[i]+y2[i]*y2[i]);
158   }
159    
160   Int_t nbinxy=400;
161   Int_t nbinz=400;
162   TH1F *hxv=new TH1F("hxv","",nbinxy,-5,5);
163   TH1F *hyv=new TH1F("hyv","",nbinxy,-5,5);
164   TH1F *hzv=new TH1F("hzv","",nbinz,-15,15);
165
166   Double_t dphi;
167   for(Int_t j=0; j<np1; j++) {
168     for(Int_t k=0; k<np2; k++) {
169       dphi=TMath::Abs(phi2[k]-phi1[j]);
170       if(dphi>180) dphi = 360-dphi;
171       if(dphi>fMaxDeltaPhi) continue; 
172       hzv->Fill(z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1));
173      }
174   }
175
176   // Fitting ...
177   Double_t max;
178   Int_t binMax;
179   Double_t maxcenter;
180   
181   max = hzv->GetMaximum();
182   binMax=hzv->GetMaximumBin();
183   maxcenter=hzv->GetBinCenter(binMax);
184   Double_t dxy=1.5;
185   
186   TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",maxcenter-dxy,maxcenter+dxy);  
187   fz->SetParameter(0,max);
188   fz->SetParameter(1,maxcenter);
189   fz->SetParameter(2,0.1);
190   fz->SetLineColor(kRed);
191   hzv->Fit("fz","RQ0");
192
193   if(TMath::Abs(x0)>0.4 || TMath::Abs(y0)>0.4) {
194     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");
195     Double_t position[3]={0,0,fz->GetParameter(1)};
196     Double_t resolution[3]={0,0,0};
197     Double_t snr[3]={0,0,0};
198     Char_t name[30];
199     sprintf(name,"Vertex");
200     fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
201     return fCurrentVertex;
202   }
203   
204   for(Int_t j=0; j<np1; j++) {
205     for(Int_t k=0; k<np2; k++) {
206       if(TMath::Abs((z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1))-fz->GetParameter(1))>fMaxDeltaZ) continue;
207       if(y2[k]==y1[j]) continue;
208       hxv->Fill(x0+(x2[k]-((x2[k]-x1[j])/(y2[k]-y1[j]))*y2[k]));
209       if(x2[k]==x1[j]) continue;
210       hyv->Fill(y0+(y2[k]-((y2[k]-y1[j])/(x2[k]-x1[j]))*x2[k]));
211     }
212   }
213     
214   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);  
215   fx->SetParameter(0,100);
216   Double_t dist=0.3;
217   Double_t xapprox=FindMaxAround(x0,hxv,dist);
218   AliDebug(1,Form("xapprox = %f",xapprox));
219   fx->SetParameter(1,xapprox);
220   Double_t difcentroid=0.07;
221   fx->SetParLimits(1,xapprox-difcentroid,xapprox+difcentroid);
222   fx->SetParameter(2,0.1);
223   fx->SetLineColor(kRed);
224   hxv->Fit("fx","RQW0"); 
225   
226   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);  
227   fy->SetParameter(0,100);
228   Double_t yapprox=FindMaxAround(y0,hyv,dist);
229   AliDebug(1,Form("yapprox = %f",yapprox));
230   fy->SetParameter(1,yapprox);
231   fy->SetParLimits(1,yapprox-difcentroid,yapprox+difcentroid);
232   fy->SetParameter(2,0.1);
233   fy->SetLineColor(kRed);
234   hyv->Fit("fy","RQW0");
235   
236   delete [] z1;
237   delete [] z2;
238   delete [] y1;
239   delete [] y2;
240   delete [] x1;
241   delete [] x2;
242   delete [] r1;
243   delete [] r2;
244   delete [] phi1;
245   delete [] phi2;
246   delete hxv;
247   delete hyv;
248   delete hzv;
249   
250   Double_t position[3]={fx->GetParameter(1),fy->GetParameter(1),fz->GetParameter(1)};
251   Double_t resolution[3]={fx->GetParameter(2),fy->GetParameter(2),fz->GetParameter(2)};
252   Double_t snr[3]={0,0,0};
253   
254   Char_t name[30];
255   sprintf(name,"Vertex");
256   fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
257
258   return fCurrentVertex;
259 }
260
261 //______________________________________________________________________
262 void AliITSVertexerIons::PhiFunc(Double_t &x,Double_t &y,Double_t &phi) {
263   // Method for the computation of the phi angle given x and y. The result is in degrees. 
264   if(y>0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578);
265   if(y>0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180;
266   if(y<0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180;
267   if(y<0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+360;;
268 }
269
270 //________________________________________________________
271 void AliITSVertexerIons::PrintStatus() const {
272   // Print current status
273   cout <<"=======================================================\n";
274   if(fCurrentVertex)fCurrentVertex->PrintStatus();
275 }
276
277 //________________________________________________________
278 Double_t AliITSVertexerIons::FindMaxAround(Double_t point, TH1F *h, Double_t distance) {
279   // It finds the maximum of h within point-distance and distance+point. 
280   Int_t maxcontent=0;
281   Int_t maxbin=0;
282   for(Int_t i=0;i<h->GetNbinsX();i++) {
283     Int_t content=(Int_t)h->GetBinContent(i);
284     Double_t center=(Double_t)h->GetBinCenter(i);
285     if(TMath::Abs(center-point)>distance) continue;
286     if(content>maxcontent) {maxcontent=content;maxbin=i;}
287   }
288   Double_t max=h->GetBinCenter(maxbin);
289   return max;
290 }
291
292