Coding rule violations corrected.
[u/mrichter/AliRoot.git] / ITS / AliITSVertexerIons.cxx
CommitLineData
b5d712e2 1#include "AliITSVertexerIons.h"
2#include "AliITSVertexerPPZ.h"
3#include "AliESDVertex.h"
c5f0f3c1 4#include "AliRun.h"
5#include "AliITS.h"
6#include "AliITSgeom.h"
88cb7938 7#include "AliITSLoader.h"
c5f0f3c1 8#include "AliITSRecPoint.h"
b5d712e2 9#include <TMath.h>
10#include <TTree.h>
c5f0f3c1 11#include <TH1.h>
12#include <TF1.h>
c5f0f3c1 13#include <TObjArray.h>
b5d712e2 14
c5f0f3c1 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 //
c5f0f3c1 24// Franco.Riggi@ct.infn.it //
25// //
fe4280bf 26// Release date: Mar 2004 //
c5f0f3c1 27// //
28// //
29//////////////////////////////////////////////////////////////////////
30
c5f0f3c1 31
b5d712e2 32ClassImp(AliITSVertexerIons)
33
fe4280bf 34//______________________________________________________________________
88cb7938 35 AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer() {
c5f0f3c1 36 // Default Constructor
37
38 fITS = 0;
83028144 39 SetNpThreshold();
fe4280bf 40 SetMaxDeltaPhi();
41 SetMaxDeltaZ();
c5f0f3c1 42}
43
44//______________________________________________________________________
88cb7938 45AliITSVertexerIons::AliITSVertexerIons(TString fn):AliITSVertexer(fn) {
c5f0f3c1 46 // Standard constructor
88cb7938 47
c5f0f3c1 48 fITS = 0;
83028144 49 SetNpThreshold();
fe4280bf 50 SetMaxDeltaPhi();
b5d712e2 51 SetMaxDeltaZ();
52}
53
54//______________________________________________________________________
55AliITSVertexerIons::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//_________________________________________________________________________
62AliITSVertexerIons& 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}
c5f0f3c1 68
69
70//______________________________________________________________________
71AliITSVertexerIons::~AliITSVertexerIons() {
72 // Default Destructor
73 fITS = 0;
74}
75
76//______________________________________________________________________
b5d712e2 77AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
78// Defines the AliESDVertex for the current event
fe4280bf 79
c5f0f3c1 80 fCurrentVertex = 0;
fe4280bf 81
88cb7938 82 AliRunLoader *rl = AliRunLoader::GetRunLoader();
83028144 83 AliITSLoader* itsloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
fe4280bf 84 if(!fITS) {
c5f0f3c1 85 fITS =(AliITS *)gAlice->GetDetector("ITS");
86 if(!fITS) {
87 Error("FindVertexForCurrentEvent","AliITS object was not found");
88 return fCurrentVertex;
89 }
90 }
88cb7938 91 fITS->SetTreeAddress();
fe4280bf 92 AliITSgeom *g2 = fITS->GetITSgeom();
88cb7938 93 TClonesArray *recpoints = fITS->RecPoints();
94 AliITSRecPoint *pnt;
fe4280bf 95 TTree *tr = itsloader->TreeR();
88cb7938 96
fe4280bf 97 Int_t npoints=0;
98 Int_t nopoints1=40000;
99 Int_t nopoints2=40000;
88cb7938 100 Float_t l[3], p[3];
c5f0f3c1 101
83028144 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];
c5f0f3c1 113
fe4280bf 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++;
2257f27e 142 }
fe4280bf 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++;
88cb7938 148 }
149 }
fe4280bf 150 }
c5f0f3c1 151
83028144 152 if(np1<fNpThreshold) {
83028144 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);
fe4280bf 155 AliITSVertexerPPZ *dovert = new AliITSVertexerPPZ("default");
83028144 156 fCurrentVertex =dovert->FindVertexForCurrentEvent(rl->GetEventNumber());
157 delete dovert;
158 return fCurrentVertex;
159 }
160
fe4280bf 161 if(!np1 || !np2) {
162 Error("FindVertexForCurrentEvent","No points in the pixer layers");
83028144 163 return fCurrentVertex;
fe4280bf 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]);
88cb7938 180 }
fe4280bf 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]);
83028144 186 }
fe4280bf 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]);
88cb7938 198 if(dphi>180) dphi = 360-dphi;
fe4280bf 199 if(dphi>fMaxDeltaPhi) continue;
200 hzv->Fill(z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1));
201 }
88cb7938 202 }
fe4280bf 203
204 // Fitting ...
205 Double_t max;
206 Int_t bin_max;
b5d712e2 207 Double_t maxcenter;
83028144 208
fe4280bf 209 max = hzv->GetMaximum();
210 bin_max=hzv->GetMaximumBin();
b5d712e2 211 maxcenter=hzv->GetBinCenter(bin_max);
fe4280bf 212 Double_t dxy=1.5;
83028144 213
b5d712e2 214 TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",maxcenter-dxy,maxcenter+dxy);
fe4280bf 215 fz->SetParameter(0,max);
b5d712e2 216 fz->SetParameter(1,maxcenter);
fe4280bf 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);
83028144 230 return fCurrentVertex;
88cb7938 231 }
83028144 232
fe4280bf 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]));
88cb7938 240 }
241 }
fe4280bf 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;
b5d712e2 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);
fe4280bf 251 fx->SetParameter(2,0.1);
252 fx->SetLineColor(kRed);
253 hxv->Fit("fx","RQW0");
83028144 254
fe4280bf 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);
b5d712e2 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);
fe4280bf 261 fy->SetParameter(2,0.1);
262 fy->SetLineColor(kRed);
263 hyv->Fit("fy","RQW0");
83028144 264
83028144 265 delete [] z1;
266 delete [] z2;
267 delete [] y1;
268 delete [] y2;
269 delete [] x1;
270 delete [] x2;
88cb7938 271 delete [] r1;
272 delete [] r2;
273 delete [] phi1;
274 delete [] phi2;
fe4280bf 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
88cb7938 283 Char_t name[30];
088e0b8d 284 if(fDebug>0)Info("FindVertexForCurrentEvent","Vertex found for event %d",evnumber);
88cb7938 285 sprintf(name,"Vertex");
d681bb2d 286 fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
88cb7938 287 return fCurrentVertex;
c5f0f3c1 288}
289
290//______________________________________________________________________
fe4280bf 291void AliITSVertexerIons::PhiFunc(Double_t &x,Double_t &y,Double_t &phi) {
b5d712e2 292 // Method for the computation of the phi angle given x and y. The result is in degrees.
fe4280bf 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;;
c5f0f3c1 297}
298
299//______________________________________________________________________
300void AliITSVertexerIons::FindVertices(){
301 // computes the vertices of the events in the range FirstEvent - LastEvent
88cb7938 302 AliRunLoader *rl = AliRunLoader::GetRunLoader();
83028144 303 AliITSLoader* itsloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
304 itsloader->LoadRecPoints("read");
ab6a6f40 305 for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
88cb7938 306 rl->GetEvent(i);
c5f0f3c1 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 }
c5f0f3c1 317}
318
319//________________________________________________________
320void AliITSVertexerIons::PrintStatus() const {
321 // Print current status
88cb7938 322 cout <<"=======================================================\n";
323 cout <<" Debug flag: "<<fDebug<<endl;
c5f0f3c1 324 cout<<"First event to be processed "<<fFirstEvent;
325 cout<<"\n Last event to be processed "<<fLastEvent<<endl;
326 if(fCurrentVertex)fCurrentVertex->PrintStatus();
327}
fe4280bf 328
329//________________________________________________________
330Double_t AliITSVertexerIons::FindMaxAround(Double_t point, TH1F *h, Double_t distance) {
b5d712e2 331 // It finds the maximum of h within point-distance and distance+point.
332 Int_t maxcontent=0;
333 Int_t maxbin=0;
fe4280bf 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);
1360bc6f 337 if(TMath::Abs(center-point)>distance) continue;
b5d712e2 338 if(content>maxcontent) {maxcontent=content;maxbin=i;}
fe4280bf 339 }
b5d712e2 340 Double_t max=h->GetBinCenter(maxbin);
fe4280bf 341 return max;
342}
343
344