]>
Commit | Line | Data |
---|---|---|
c5f0f3c1 | 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" | |
88cb7938 | 14 | #include "AliITSLoader.h" |
c5f0f3c1 | 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> | |
d681bb2d | 22 | #include <AliESDVertex.h> |
c5f0f3c1 | 23 | #include <TObjArray.h> |
24 | #include <TObject.h> | |
83028144 | 25 | #include <AliITSVertexerPPZ.h> |
c5f0f3c1 | 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 // | |
83028144 | 34 | // |
c5f0f3c1 | 35 | // Giuseppe.Lore@ct.infn.it // |
c5f0f3c1 | 36 | // Franco.Riggi@ct.infn.it // |
37 | // // | |
38 | // Release date: May 2001 // | |
39 | // // | |
40 | // // | |
41 | ////////////////////////////////////////////////////////////////////// | |
42 | ||
43 | ClassImp(AliITSVertexerIons) | |
44 | ||
45 | ||
46 | ||
83028144 | 47 | //______________________________________________________________________ |
88cb7938 | 48 | AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer() { |
c5f0f3c1 | 49 | // Default Constructor |
50 | ||
51 | fITS = 0; | |
83028144 | 52 | SetNpThreshold(); |
c5f0f3c1 | 53 | } |
54 | ||
55 | //______________________________________________________________________ | |
88cb7938 | 56 | AliITSVertexerIons::AliITSVertexerIons(TString fn):AliITSVertexer(fn) { |
c5f0f3c1 | 57 | // Standard constructor |
88cb7938 | 58 | |
c5f0f3c1 | 59 | fITS = 0; |
83028144 | 60 | SetNpThreshold(); |
c5f0f3c1 | 61 | } |
62 | ||
63 | ||
64 | //______________________________________________________________________ | |
65 | AliITSVertexerIons::~AliITSVertexerIons() { | |
66 | // Default Destructor | |
67 | fITS = 0; | |
68 | } | |
69 | ||
70 | //______________________________________________________________________ | |
d681bb2d | 71 | AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){ |
72 | // Defines the AliESDVertex for the current event | |
c5f0f3c1 | 73 | fCurrentVertex = 0; |
83028144 | 74 | Double_t position[3]; |
75 | Double_t resolution[3]; | |
76 | Double_t snr[3]; | |
88cb7938 | 77 | AliRunLoader *rl = AliRunLoader::GetRunLoader(); |
83028144 | 78 | AliITSLoader* itsloader = (AliITSLoader*) rl->GetLoader("ITSLoader"); |
c5f0f3c1 | 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 | } | |
88cb7938 | 86 | fITS->SetTreeAddress(); |
87 | AliITSgeom *g2 = fITS->GetITSgeom(); | |
88 | TClonesArray *recpoints = fITS->RecPoints(); | |
89 | AliITSRecPoint *pnt; | |
90 | ||
83028144 | 91 | Int_t nopoints1 = 0; |
92 | Int_t nopoints2 = 0; | |
93 | Double_t vzero[3]; | |
94 | Double_t aspar[2]; | |
c5f0f3c1 | 95 | |
88cb7938 | 96 | //------------ Rough Vertex evaluation --------------------------------- |
c5f0f3c1 | 97 | |
83028144 | 98 | Int_t i,npoints,ipoint,j,k,max,binmax; |
99 | Double_t zCentroid; | |
88cb7938 | 100 | Float_t l[3], p[3]; |
c5f0f3c1 | 101 | |
88cb7938 | 102 | Double_t mxpiu = 0; |
103 | Double_t mxmeno = 0; | |
104 | Double_t mypiu = 0; | |
105 | Double_t mymeno = 0; | |
c5f0f3c1 | 106 | |
88cb7938 | 107 | TH1F *hITSz1 = new TH1F("hITSz1","",100,-14.35,14.35); |
83028144 | 108 | TTree *tr = itsloader->TreeR(); |
88cb7938 | 109 | for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++) |
110 | { | |
111 | fITS->ResetRecPoints(); | |
83028144 | 112 | tr->GetEvent(i); |
88cb7938 | 113 | npoints = recpoints->GetEntries(); |
114 | for (ipoint=0;ipoint<npoints;ipoint++) { | |
115 | pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint); | |
116 | l[0]=pnt->GetX(); | |
117 | l[1]=0; | |
118 | l[2]=pnt->GetZ(); | |
119 | g2->LtoG(i, l, p); | |
120 | if(i<80 && TMath::Abs(p[2])<14.35) { | |
121 | if(p[0]>0) mxpiu++; | |
122 | if(p[0]<0) mxmeno++; | |
123 | if(p[1]>0) mypiu++; | |
124 | if(p[1]<0) mymeno++; | |
125 | hITSz1->Fill(p[2]); | |
126 | } | |
83028144 | 127 | if(i>=80 && TMath::Abs(p[2])<14.35) nopoints2++; |
88cb7938 | 128 | } |
129 | } | |
c5f0f3c1 | 130 | |
83028144 | 131 | nopoints1 = (Int_t)(hITSz1->GetEntries()); |
c5f0f3c1 | 132 | |
83028144 | 133 | aspar[0] = (mxpiu-mxmeno)/(mxpiu+mxmeno); |
134 | aspar[1] = (mypiu-mymeno)/(mypiu+mymeno); | |
c5f0f3c1 | 135 | |
83028144 | 136 | vzero[0] = 5.24441*aspar[0]; |
137 | vzero[1] = 5.24441*aspar[1]; | |
c5f0f3c1 | 138 | |
83028144 | 139 | zCentroid= TMath::Abs(hITSz1->GetMean()); |
140 | vzero[2] = -5.31040e-02+1.42198*zCentroid+7.44718e-01*TMath::Power(zCentroid,2) | |
141 | -5.73426e-01*TMath::Power(zCentroid,3)+2.01500e-01*TMath::Power(zCentroid,4) | |
142 | -3.34118e-02*TMath::Power(zCentroid,5)+2.20816e-03*TMath::Power(zCentroid,6); | |
c5f0f3c1 | 143 | |
83028144 | 144 | if(hITSz1->GetMean()<0) vzero[2] = -vzero[2]; |
c5f0f3c1 | 145 | |
e83fe598 | 146 | //cout << "\nXvzero: " << vzero[0] << " cm" << ""; |
147 | //cout << "\nYvzero: " << vzero[1] << " cm" << ""; | |
148 | //cout << "\nZvzero: " << vzero[2] << " cm" << "\n"; | |
c5f0f3c1 | 149 | |
88cb7938 | 150 | delete hITSz1; |
c5f0f3c1 | 151 | |
83028144 | 152 | Double_t dphi,r,deltaZ,a,b; |
88cb7938 | 153 | Int_t np1=0; |
154 | Int_t np2=0; | |
155 | Int_t niter=0; | |
c5f0f3c1 | 156 | |
83028144 | 157 | Double_t deltaPhiZ = 0.08; |
c5f0f3c1 | 158 | |
83028144 | 159 | Float_t mag[3]; |
88cb7938 | 160 | Float_t origin[3]; |
161 | for(Int_t okm=0;okm<3;okm++) origin[okm]=0; | |
83028144 | 162 | gAlice->Field()->Field(origin,mag); |
c5f0f3c1 | 163 | |
e83fe598 | 164 | deltaPhiZ = deltaPhiZ*TMath::Abs(mag[2])/2; |
83028144 | 165 | Double_t deltaPhiXY = 1.0; |
c5f0f3c1 | 166 | |
e83fe598 | 167 | //cout << "\ndeltaPhiZ: " << deltaPhiZ << " deg" << "\n"; |
168 | //cout << "deltaPhiXY: " << deltaPhiXY << " deg" << "\n"; | |
c5f0f3c1 | 169 | |
83028144 | 170 | Double_t *z1, *z2, *y1, *y2, *x1, *x2, *phi1, *phi2, *r1, *r2; |
171 | z1=new Double_t[nopoints1]; | |
172 | z2=new Double_t[nopoints2]; | |
173 | y1=new Double_t[nopoints1]; | |
174 | y2=new Double_t[nopoints2]; | |
175 | x1=new Double_t[nopoints1]; | |
176 | x2=new Double_t[nopoints2]; | |
177 | phi1=new Double_t[nopoints1]; | |
178 | phi2=new Double_t[nopoints2]; | |
179 | r1=new Double_t[nopoints1]; | |
180 | r2=new Double_t[nopoints2]; | |
c5f0f3c1 | 181 | |
83028144 | 182 | deltaZ = 4.91617e-01+2.67567e-02*vzero[2]+1.49626e-02*TMath::Power(vzero[2],2); |
183 | Float_t multFactorZ = 28000./(Float_t)nopoints1; | |
184 | Int_t nbin=(Int_t)((deltaZ/0.005)/multFactorZ); | |
88cb7938 | 185 | Int_t nbinxy=250; |
83028144 | 186 | Int_t *vectorBinZ,*vectorBinXY; |
187 | vectorBinZ=new Int_t[nbin]; | |
188 | vectorBinXY=new Int_t[nbinxy]; | |
88cb7938 | 189 | Float_t f1= 0; |
190 | Float_t f2= 0; | |
83028144 | 191 | Double_t sigma,averagebg; |
c5f0f3c1 | 192 | |
83028144 | 193 | TH1D *hITSZv = new TH1D("hITSZv","",nbin,vzero[2]-deltaZ,vzero[2]+deltaZ); |
88cb7938 | 194 | TH1D *hITSXv = new TH1D("hITSXv","",nbinxy,-3,3); |
195 | TH1D *hITSYv = new TH1D("hITSYv","",nbinxy,-3,3); | |
c5f0f3c1 | 196 | |
e83fe598 | 197 | //cout << "deltaZeta: " << deltaZ << " cm" << "\n"; |
c5f0f3c1 | 198 | |
199 | ||
88cb7938 | 200 | start: |
c5f0f3c1 | 201 | |
83028144 | 202 | hITSZv->Add(hITSZv,-1); |
203 | hITSXv->Add(hITSXv,-1); | |
204 | hITSYv->Add(hITSYv,-1); | |
c5f0f3c1 | 205 | |
88cb7938 | 206 | np1=np2=0; |
c5f0f3c1 | 207 | |
c5f0f3c1 | 208 | |
88cb7938 | 209 | for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++) |
210 | { | |
211 | fITS->ResetRecPoints(); | |
83028144 | 212 | itsloader->TreeR()->GetEvent(i); |
88cb7938 | 213 | npoints = recpoints->GetEntries(); |
214 | for (ipoint=0;ipoint<npoints;ipoint++) { | |
c5f0f3c1 | 215 | |
88cb7938 | 216 | pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint); |
217 | l[0]=pnt->GetX(); | |
218 | l[1]=0; | |
219 | l[2]=pnt->GetZ(); | |
220 | g2->LtoG(i, l, p); | |
c5f0f3c1 | 221 | |
88cb7938 | 222 | if(i<80 && TMath::Abs(p[2])<14.35) { |
83028144 | 223 | p[0]=p[0]-vzero[0]; |
224 | p[1]=p[1]-vzero[1]; | |
88cb7938 | 225 | r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2)); |
83028144 | 226 | y1[np1]=p[1]; |
227 | x1[np1]=p[0]; | |
228 | z1[np1]=p[2]; | |
88cb7938 | 229 | r1[np1]=r; |
230 | phi1[np1]=PhiFunc(p); | |
231 | np1++; | |
232 | } | |
233 | ||
234 | if(i>=80 && TMath::Abs(p[2])<14.35) { | |
83028144 | 235 | p[0]=p[0]-vzero[0]; |
236 | p[1]=p[1]-vzero[1]; | |
88cb7938 | 237 | r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2)); |
83028144 | 238 | y2[np2]=p[1]; |
239 | x2[np2]=p[0]; | |
240 | z2[np2]=p[2]; | |
88cb7938 | 241 | r2[np2]=r; |
242 | phi2[np2]=PhiFunc(p); | |
243 | np2++; | |
244 | } | |
c5f0f3c1 | 245 | |
88cb7938 | 246 | } |
247 | } | |
c5f0f3c1 | 248 | |
88cb7938 | 249 | //------------------ Correlation between rec points ---------------------- |
c5f0f3c1 | 250 | |
83028144 | 251 | |
252 | if(np1<fNpThreshold) { | |
253 | /* cout << "Warning: AliITSVertexerIons finder is not reliable for low multiplicity events. May be you have to use AliITSVertexerPPZ.\n"; | |
254 | position[0]=vzero[0]; | |
255 | position[1]=vzero[1]; | |
256 | position[2]=vzero[2]; | |
257 | resolution[0]=-123; | |
258 | resolution[1]=-123; | |
259 | resolution[2]=-123; | |
260 | snr[0]=-123; | |
261 | snr[1]=-123; | |
262 | snr[2]=-123; | |
263 | Char_t name[30]; | |
264 | sprintf(name,"Vertex"); | |
d681bb2d | 265 | fCurrentVertex = new AliESDVertex(position,resolution,snr,name); |
83028144 | 266 | return fCurrentVertex;*/ |
267 | ||
268 | Warning("FindVertexForCurrentEvent","AliITSVertexerIons finder is not reliable for low multiplicity events. Switching to AliITSVertexerPPZ with default parameters...\n"); | |
269 | Warning("FindVertexForCurrentEvent","N rec points = %d - Threshold is %d",np1,fNpThreshold); | |
270 | AliITSVertexerPPZ *dovert = new AliITSVertexerPPZ("default"); | |
271 | fCurrentVertex =dovert->FindVertexForCurrentEvent(rl->GetEventNumber()); | |
272 | delete dovert; | |
273 | return fCurrentVertex; | |
274 | } | |
275 | ||
88cb7938 | 276 | |
277 | for(j=0; j<(np2)-1; j++) { | |
278 | for(k=0; k<(np1)-1; k++) { | |
279 | dphi=TMath::Abs(phi1[k]-phi2[j]); | |
280 | if(dphi>180) dphi = 360-dphi; | |
83028144 | 281 | if(dphi<deltaPhiZ && TMath::Abs((z1[k]-(z2[j]-z1[k])/((r2[j]/r1[k])-1))-vzero[2]) |
282 | <deltaZ) hITSZv->Fill(z1[k]-(z2[j]-z1[k])/((r2[j]/r1[k])-1)); | |
88cb7938 | 283 | } |
284 | } | |
c5f0f3c1 | 285 | |
88cb7938 | 286 | //cout << "\nNumber of used pairs: \n"; |
83028144 | 287 | |
288 | a = vzero[2]-deltaZ; | |
289 | b = vzero[2]+deltaZ; | |
88cb7938 | 290 | max=(Int_t) hITSZv->GetMaximum(); |
83028144 | 291 | binmax=hITSZv->GetMaximumBin(); |
88cb7938 | 292 | sigma=0; |
83028144 | 293 | for(i=0;i<nbin;i++) vectorBinZ[i]=(Int_t)hITSZv->GetBinContent(i); |
294 | for(i=0;i<10;i++) f1=f1+vectorBinZ[i]/10; | |
295 | for(i=nbin-10;i<nbin;i++) f2=f2+vectorBinZ[i]/10; | |
296 | averagebg=(f1+f2)/2; | |
88cb7938 | 297 | for(i=0;i<nbin;i++) { |
83028144 | 298 | if(vectorBinZ[i]-averagebg>(max-averagebg)*0.4 && |
299 | vectorBinZ[i]-averagebg<(max-averagebg)*0.7) { | |
300 | sigma=hITSZv->GetBinCenter(binmax)-hITSZv->GetBinCenter(i); | |
88cb7938 | 301 | sigma=TMath::Abs(sigma); |
302 | if(sigma==0) sigma=0.05; | |
303 | } | |
304 | } | |
83028144 | 305 | |
306 | ||
88cb7938 | 307 | TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",a,b); |
83028144 | 308 | |
88cb7938 | 309 | fz->SetParameter(0,max); |
83028144 | 310 | if(niter==0) {Double_t temp = hITSZv->GetBinCenter(binmax); vzero[2]=temp;} |
311 | fz->SetParameter(1,vzero[2]); | |
88cb7938 | 312 | fz->SetParameter(2,sigma); |
83028144 | 313 | fz->SetParameter(3,averagebg); |
88cb7938 | 314 | fz->SetParLimits(2,0,999); |
315 | fz->SetParLimits(3,0,999); | |
83028144 | 316 | |
317 | hITSZv->Fit("fz","RQ0"); | |
318 | ||
319 | snr[2] = fz->GetParameter(0)/fz->GetParameter(3); | |
320 | if(snr[2]<0.) { | |
321 | Error("FindVertexForCurrentEvent","\nNegative Signal to noise ratio for z!!!\n"); | |
322 | Error("FindVertexForCurrentEvent","The algorithm cannot find the z vertex position.\n"); | |
323 | // exit(123456789); | |
324 | return fCurrentVertex; | |
88cb7938 | 325 | } |
83028144 | 326 | else { |
327 | position[2] = fz->GetParameter(1); | |
328 | if(position[2]<0) | |
329 | { | |
330 | position[2]=position[2]-TMath::Abs(position[2])*1.11/10000; | |
331 | } | |
332 | else | |
333 | { | |
334 | position[2]=position[2]+TMath::Abs(position[2])*1.11/10000; | |
335 | } | |
336 | } | |
337 | resolution[2] = fz->GetParError(1); | |
338 | ||
88cb7938 | 339 | for(j=0; j<(np2)-1; j++) { |
340 | for(k=0; k<(np1)-1; k++) { | |
341 | dphi=TMath::Abs(phi1[k]-phi2[j]); | |
342 | if(dphi>180) dphi = 360-dphi; | |
83028144 | 343 | if(dphi>deltaPhiXY) continue; |
344 | if(TMath::Abs((z1[k]-(z2[j]-z1[k])/((r2[j]/r1[k])-1))-position[2]) | |
345 | <4*resolution[2]) { | |
346 | hITSXv->Fill(vzero[0]+(x2[j]-((x2[j]-x1[k])/(y2[j]-y1[k]))*y2[j])); | |
347 | hITSYv->Fill(vzero[1]+(y2[j]-((y2[j]-y1[k])/(x2[j]-x1[k]))*x2[j])); | |
c5f0f3c1 | 348 | } |
88cb7938 | 349 | } |
350 | } | |
83028144 | 351 | |
352 | TF1 *fx = new TF1("fx","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",vzero[0]-0.5,vzero[0]+0.5); | |
353 | ||
88cb7938 | 354 | max=(Int_t) hITSXv->GetMaximum(); |
83028144 | 355 | binmax=hITSXv->GetMaximumBin(); |
88cb7938 | 356 | sigma=0; |
357 | f1=f2=0; | |
83028144 | 358 | for(i=0;i<nbinxy;i++) vectorBinXY[i]=(Int_t)hITSXv->GetBinContent(i); |
359 | for(i=0;i<10;i++) f1=f1+vectorBinXY[i]/10; | |
360 | for(i=nbinxy-10;i<nbinxy;i++) f2=f2+vectorBinXY[i]/10; | |
361 | averagebg=(f1+f2)/2; | |
88cb7938 | 362 | for(i=0;i<nbinxy;i++) { |
83028144 | 363 | if(vectorBinXY[i]-averagebg>(max-averagebg)*0.4 && |
364 | vectorBinXY[i]-averagebg<(max-averagebg)*0.7) { | |
365 | sigma=hITSXv->GetBinCenter(binmax)-hITSXv->GetBinCenter(i); | |
88cb7938 | 366 | sigma=TMath::Abs(sigma); |
367 | if(sigma==0) sigma=0.05; | |
368 | } | |
369 | } | |
83028144 | 370 | |
371 | ||
88cb7938 | 372 | fx->SetParameter(0,max); |
83028144 | 373 | fx->SetParameter(1,vzero[0]); |
88cb7938 | 374 | fx->SetParameter(2,sigma); |
83028144 | 375 | fx->SetParameter(3,averagebg); |
376 | ||
377 | hITSXv->Fit("fx","RQ0"); | |
378 | ||
379 | snr[0] = fx->GetParameter(0)/fx->GetParameter(3); | |
380 | ||
381 | if(snr[0]<0.) { | |
382 | Error("FindVertexForCurrentEvent","\nNegative Signal to noise ratio for x!\n"); | |
383 | Error("FindVertexForCurrentEvent","The algorithm cannot find the x vertex position\n"); | |
384 | // exit(123456789); | |
385 | return fCurrentVertex; | |
88cb7938 | 386 | } |
83028144 | 387 | else { |
388 | position[0]=fx->GetParameter(1); | |
389 | resolution[0]=fx->GetParError(1); | |
390 | } | |
391 | ||
392 | TF1 *fy = new TF1("fy","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",vzero[1]-0.5,vzero[1]+0.5); | |
393 | ||
88cb7938 | 394 | max=(Int_t) hITSYv->GetMaximum(); |
83028144 | 395 | binmax=hITSYv->GetMaximumBin(); |
88cb7938 | 396 | sigma=0; |
397 | f1=f2=0; | |
83028144 | 398 | for(i=0;i<nbinxy;i++) vectorBinXY[i]=(Int_t)hITSYv->GetBinContent(i); |
399 | for(i=0;i<10;i++) f1=f1+vectorBinXY[i]/10; | |
400 | for(i=nbinxy-10;i<nbinxy;i++) f2=f2+vectorBinXY[i]/10; | |
401 | averagebg=(f1+f2)/2; | |
88cb7938 | 402 | for(i=0;i<nbinxy;i++) { |
83028144 | 403 | if(vectorBinXY[i]-averagebg>(max-averagebg)*0.4 && |
404 | vectorBinXY[i]-averagebg<(max-averagebg)*0.7) { | |
405 | sigma=hITSYv->GetBinCenter(binmax)-hITSYv->GetBinCenter(i); | |
88cb7938 | 406 | sigma=TMath::Abs(sigma); |
407 | if(sigma==0) sigma=0.05; | |
408 | } | |
409 | } | |
83028144 | 410 | |
88cb7938 | 411 | fy->SetParameter(0,max); |
83028144 | 412 | fy->SetParameter(1,vzero[1]); |
88cb7938 | 413 | fy->SetParameter(2,sigma); |
83028144 | 414 | fy->SetParameter(3,averagebg); |
415 | ||
416 | hITSYv->Fit("fy","RQ0"); | |
417 | ||
418 | snr[1] = fy->GetParameter(0)/fy->GetParameter(3); | |
419 | if(snr[1]<0.) { | |
420 | Error("FindVertexForCurrentEvent","\nNegative Signal to noise ratio for y!\n"); | |
421 | Error("FindVertexForCurrentEvent","The algorithm cannot find the y vertex position.\n"); | |
422 | // exit(123456789); | |
423 | return fCurrentVertex; | |
88cb7938 | 424 | } |
83028144 | 425 | else { |
426 | position[1]=fy->GetParameter(1); | |
427 | resolution[1]=fy->GetParError(1); | |
428 | } | |
429 | ||
430 | ||
431 | /*cout << "iter = " << niter << " -> x = " << position[0] << endl; | |
432 | cout << "iter = " << niter << " -> y = " << position[1] << endl; | |
433 | cout << "iter = " << niter << " -> z = " << position[2] << endl; | |
434 | cout << "---\n";*/ | |
88cb7938 | 435 | |
436 | niter++; | |
c5f0f3c1 | 437 | |
83028144 | 438 | vzero[0] = position[0]; |
439 | vzero[1] = position[1]; | |
440 | vzero[2] = position[2]; | |
c5f0f3c1 | 441 | |
88cb7938 | 442 | if(niter<3) goto start; // number of iterations |
c5f0f3c1 | 443 | |
444 | ||
88cb7938 | 445 | cout<<"Number of iterations: "<<niter<<endl; |
83028144 | 446 | cout << "\nNo. of Points on the two pixel layers: "<<np1<<" "<<np2<<endl; |
c5f0f3c1 | 447 | |
83028144 | 448 | delete [] z1; |
449 | delete [] z2; | |
450 | delete [] y1; | |
451 | delete [] y2; | |
452 | delete [] x1; | |
453 | delete [] x2; | |
88cb7938 | 454 | delete [] r1; |
455 | delete [] r2; | |
456 | delete [] phi1; | |
457 | delete [] phi2; | |
83028144 | 458 | delete [] vectorBinZ; |
459 | delete [] vectorBinXY; | |
c5f0f3c1 | 460 | |
88cb7938 | 461 | delete hITSZv; |
462 | delete hITSXv; | |
463 | delete hITSYv; | |
464 | Char_t name[30]; | |
465 | // sprintf(name,"Vertex_%d",evnumber); | |
088e0b8d | 466 | if(fDebug>0)Info("FindVertexForCurrentEvent","Vertex found for event %d",evnumber); |
88cb7938 | 467 | sprintf(name,"Vertex"); |
d681bb2d | 468 | fCurrentVertex = new AliESDVertex(position,resolution,snr,name); |
88cb7938 | 469 | return fCurrentVertex; |
c5f0f3c1 | 470 | } |
471 | ||
472 | //______________________________________________________________________ | |
473 | Double_t AliITSVertexerIons::PhiFunc(Float_t p[]) { | |
88cb7938 | 474 | Double_t phi=0; |
c5f0f3c1 | 475 | |
88cb7938 | 476 | if(p[1]>0 && p[0]>0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578); |
477 | if(p[1]>0 && p[0]<0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+180; | |
478 | if(p[1]<0 && p[0]<0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+180; | |
479 | if(p[1]<0 && p[0]>0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+360; | |
480 | return phi; | |
c5f0f3c1 | 481 | } |
482 | ||
483 | //______________________________________________________________________ | |
484 | void AliITSVertexerIons::FindVertices(){ | |
485 | // computes the vertices of the events in the range FirstEvent - LastEvent | |
88cb7938 | 486 | AliRunLoader *rl = AliRunLoader::GetRunLoader(); |
83028144 | 487 | AliITSLoader* itsloader = (AliITSLoader*) rl->GetLoader("ITSLoader"); |
488 | itsloader->LoadRecPoints("read"); | |
ab6a6f40 | 489 | for(Int_t i=fFirstEvent;i<=fLastEvent;i++){ |
88cb7938 | 490 | rl->GetEvent(i); |
c5f0f3c1 | 491 | FindVertexForCurrentEvent(i); |
492 | if(fCurrentVertex){ | |
493 | WriteCurrentVertex(); | |
494 | } | |
495 | else { | |
496 | if(fDebug>0){ | |
497 | cout<<"Vertex not found for event "<<i<<endl; | |
498 | } | |
499 | } | |
500 | } | |
c5f0f3c1 | 501 | } |
502 | ||
503 | //________________________________________________________ | |
504 | void AliITSVertexerIons::PrintStatus() const { | |
505 | // Print current status | |
88cb7938 | 506 | cout <<"=======================================================\n"; |
507 | cout <<" Debug flag: "<<fDebug<<endl; | |
c5f0f3c1 | 508 | cout<<"First event to be processed "<<fFirstEvent; |
509 | cout<<"\n Last event to be processed "<<fLastEvent<<endl; | |
510 | if(fCurrentVertex)fCurrentVertex->PrintStatus(); | |
511 | } |