]>
Commit | Line | Data |
---|---|---|
eb35e591 | 1 | /************************************************************************** |
2 | * Copyright(c) 2006-2008, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | //----------------------------------------------------------------- | |
17 | // AliITSVertexer3DTapan class | |
18 | // This is a class for the 3d vertex finding | |
19 | // Origin: Tapan Nayak, VECC-CERN, Tapan.Nayak@cern.ch | |
20 | //----------------------------------------------------------------- | |
21 | ||
22 | #include <TH1.h> | |
23 | #include <TTree.h> | |
24 | #include <TClonesArray.h> | |
25 | ||
26 | #include <AliITSVertexer3DTapan.h> | |
27 | #include <AliITSRecPoint.h> | |
28 | #include <AliITSgeom.h> | |
29 | #include <AliESDVertex.h> | |
30 | ||
31 | ClassImp(AliITSVertexer3DTapan) | |
32 | ||
308c2f7c | 33 | void AliITSVertexer3DTapan::LoadClusters(TTree *cTree) { |
eb35e591 | 34 | //-------------------------------------------------------------------- |
35 | //This function loads the SPD clusters | |
36 | //-------------------------------------------------------------------- | |
37 | TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy; | |
38 | TBranch *branch=cTree->GetBranch("ITSRecPoints"); | |
39 | branch->SetAddress(&clusters); | |
40 | ||
41 | Int_t nentr=cTree->GetEntries(),nc1=0,nc2=0; | |
42 | for (Int_t i=0; i<nentr; i++) { | |
43 | if (!cTree->GetEvent(i)) continue; | |
44 | // | |
45 | // Below: | |
46 | // "alpha" is the angle from the global X-axis to the | |
47 | // local GEANT X'-axis ( rot[0]=cos(alpha) and rot[1]=sin(alpha) ) | |
48 | // "phi" is the angle from the global X-axis to the | |
49 | // local cluster X"-axis | |
50 | // | |
51 | ||
52 | // Double_t rot[9]; fITSgeom->GetRotMatrix(i,rot); | |
53 | Int_t lay,lad,det; fITSgeom->GetModuleId(i,lay,lad,det); | |
54 | ||
55 | if (lay>2) break; //load the SPD clusters only | |
56 | ||
57 | /* | |
58 | Float_t tx,ty,tz; fITSgeom->GetTrans(lay,lad,det,tx,ty,tz); | |
59 | ||
60 | Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi(); | |
61 | Double_t phi=TMath::Pi()/2+alpha; | |
62 | ||
63 | if (lay==1) phi+=TMath::Pi(); | |
64 | Double_t cp=TMath::Cos(phi), sp=TMath::Sin(phi); | |
65 | Double_t r=tx*cp+ty*sp; | |
66 | */ | |
67 | ||
68 | Int_t ncl=clusters->GetEntriesFast(); | |
69 | Float_t hPhi; | |
70 | while (ncl--) { | |
71 | AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl); | |
72 | Float_t pos[3]; | |
73 | c->GetGlobalXYZ(pos); | |
74 | if (lay==1) { | |
75 | /* fX1[nc1]= r*cp - c->GetY()*sp; | |
76 | fY1[nc1]= r*sp + c->GetY()*cp; | |
77 | fZ1[nc1]= c->GetZ(); */ | |
78 | fX1[nc1] = pos[0]; fY1[nc1] = pos[1]; fZ1[nc1] = pos[2]; | |
79 | CalculatePhi(fX1[nc1], fY1[nc1], hPhi); | |
80 | fPhi1[nc1]= hPhi; | |
81 | nc1++; | |
82 | } else { | |
83 | /* fX2[nc2]= r*cp - c->GetY()*sp; | |
84 | fY2[nc2]= r*sp + c->GetY()*cp; | |
85 | fZ2[nc2]= c->GetZ(); */ | |
86 | fX2[nc2] = pos[0]; fY2[nc2] = pos[1]; fZ2[nc2] = pos[2]; | |
87 | CalculatePhi(fX2[nc2], fY2[nc2], hPhi); | |
88 | fPhi2[nc2]= hPhi; | |
89 | nc2++; | |
90 | } | |
91 | } | |
92 | } | |
93 | ficlu1 = nc1; ficlu2 = nc2; | |
94 | AliInfo(Form("Number of clusters: %d (first layer) and %d (second layer)",ficlu1,ficlu2)); | |
eb35e591 | 95 | } |
96 | ||
308c2f7c | 97 | AliESDVertex *AliITSVertexer3DTapan::FindVertexForCurrentEvent(TTree *cTree) { |
eb35e591 | 98 | // |
99 | // This function reconstructs .... | |
100 | // | |
101 | // | |
308c2f7c | 102 | LoadClusters(cTree); |
eb35e591 | 103 | |
104 | Double_t pos[3], postemp[3], sigpos[3]; | |
105 | Int_t ncontr, ncontrtemp; | |
106 | Float_t cuts[3]; | |
107 | Int_t vtxstatus=0; | |
108 | ||
109 | //.... | |
110 | pos[0] = 0.; pos[1] = 0.; pos[2] = 0.; | |
111 | cuts[0]=1.; cuts[1]=1.; cuts[2]=20.; | |
112 | CalculateVertex3d1(pos, cuts, ncontr); | |
113 | if(ncontr==0) { | |
114 | pos[0] = 9999.; pos[1] = 9999.; pos[2] = 9999.; | |
115 | vtxstatus = -1; | |
116 | } | |
117 | AliInfo(Form("1st step: %d %f %f %f st=%d",ncontr,pos[0],pos[1],pos[2],vtxstatus)); | |
118 | ||
119 | if(vtxstatus == 0) { | |
120 | ncontrtemp = ncontr; postemp[0] = pos[0]; postemp[1] = pos[1]; postemp[2] = pos[2]; | |
121 | cuts[0]=0.3; cuts[1]=0.3; cuts[2]=1.; | |
122 | CalculateVertex3d1(pos, cuts, ncontr); | |
123 | if(ncontr==0) { | |
124 | ncontr = ncontrtemp; pos[0] = postemp[0]; pos[1] = postemp[1]; pos[2] = postemp[2]; | |
125 | vtxstatus = 2; | |
126 | } | |
127 | AliInfo(Form("2nd step: %d %f %f %f st=%d",ncontr,pos[0],pos[1],pos[2],vtxstatus)); | |
128 | } | |
129 | ||
130 | if(vtxstatus == 0) { | |
131 | ncontrtemp = ncontr; postemp[0] = pos[0]; postemp[1] = pos[1]; postemp[2] = pos[2]; | |
132 | cuts[0]=0.25; cuts[1]=0.25; cuts[2]=1.0; | |
133 | CalculateVertex3d2(pos, cuts, ncontr, sigpos); | |
134 | if(ncontr==0) { | |
135 | ncontr = ncontrtemp; pos[0] = postemp[0]; pos[1] = postemp[1]; pos[2] = postemp[2]; | |
136 | vtxstatus = 3; | |
137 | } | |
138 | AliInfo(Form("3rd step: %d %f %f %f st=%d",ncontr,pos[0],pos[1],pos[2],vtxstatus)); | |
139 | } | |
140 | ||
141 | if(vtxstatus == 0) { | |
142 | ncontrtemp = ncontr; postemp[0] = pos[0]; postemp[1] = pos[1]; postemp[2] = pos[2]; | |
143 | cuts[0]=0.1; cuts[1]=0.1; cuts[2]=0.2; | |
144 | CalculateVertex3d2(pos, cuts, ncontr, sigpos); | |
145 | if(ncontr==0) { | |
146 | ncontr = ncontrtemp; pos[0] = postemp[0]; pos[1] = postemp[1]; pos[2] = postemp[2]; | |
147 | vtxstatus = 4; | |
148 | } | |
149 | AliInfo(Form("4th step: %d %f %f %f st=%d",ncontr,pos[0],pos[1],pos[2],vtxstatus)); | |
150 | } | |
151 | AliInfo(Form("Final step: %d %f %f %f st=%d",ncontr,pos[0],pos[1],pos[2],vtxstatus)); | |
152 | ||
308c2f7c | 153 | return new AliESDVertex(pos,sigpos,(Double_t)vtxstatus,ncontr,"AliITSVertexer3DTapan"); |
154 | ||
eb35e591 | 155 | } |
156 | ||
157 | ||
158 | void AliITSVertexer3DTapan::CalculateVertex3d1(Double_t pos[3], Float_t cuts[3], Int_t &ncontr) { | |
159 | // | |
160 | // This function reconstructs first two steps of vertex | |
161 | // | |
162 | ||
163 | Double_t p1[4], p2[4], p3[4], p4[4]; | |
164 | Double_t pa[3], pb[3]; | |
165 | Double_t hphi1, hphi2, hphi3, hphi4; | |
166 | ||
167 | ncontr = 0; | |
168 | Float_t phicut = 1.0; | |
169 | Double_t distance; Float_t distancecut = 1.0; | |
170 | Int_t ibin=20; Float_t ilow=-1.; Float_t ihigh=1.; | |
171 | Int_t ibinz=400; Float_t ilowz=-20.; Float_t ihighz=20.; | |
172 | ||
173 | TH1F *hx = new TH1F("hx","", ibin, ilow, ihigh); | |
174 | TH1F *hy = new TH1F("hy","", ibin, ilow, ihigh); | |
175 | TH1F *hz = new TH1F("hz","", ibinz,ilowz,ihighz); | |
176 | ||
177 | for (Int_t ip1=0; ip1<ficlu1; ip1++) { | |
178 | // Two points on layer1: p1 and p3 | |
179 | p1[0] = fX1[ip1]; p1[1] = fY1[ip1]; p1[2] = fZ1[ip1]; | |
180 | p3[0] = fX1[ip1+1]; p3[1] = fY1[ip1+1]; p3[2] = fZ1[ip1+1]; | |
181 | hphi1 = fPhi1[ip1]; hphi3 = fPhi1[ip1+1]; | |
182 | ||
183 | for (Int_t ip2=0; ip2<ficlu2; ip2++) { | |
184 | // Two points on layer 2: p2 and p4 | |
185 | p2[0] = fX2[ip2]; p2[1] = fY2[ip2]; p2[2] = fZ2[ip2]; | |
186 | p4[0] = fX2[ip2+1]; p4[1] = fY2[ip2+1]; p4[2] = fZ2[ip2+1]; | |
187 | hphi2 = fPhi2[ip2]; hphi4 = fPhi2[ip2+1]; | |
188 | ||
189 | // First line is formed by p1-p2 and second line by p3-p4 | |
190 | // We find two points on each line which form the closest distance of the two lines | |
191 | // pa[0],pa[1],pa[2]: points on line 1 and pb[0],pb[1],pb[2]: points on line 2 | |
192 | // Next: Consider x, y and z to be less than cuts[0], cuts[1] and cuts[2], respectively | |
193 | ||
194 | if(TMath::Abs(hphi1-hphi2)<phicut && TMath::Abs(hphi3-hphi4)<phicut){ | |
195 | CalculateLine(p1, p2, p3, p4, pa, pb); | |
196 | ||
197 | if (pa[0]>pos[0]-cuts[0] && pa[0]<pos[0]+cuts[0] && pa[1]>pos[1]-cuts[1] && pa[1]<pos[1]+cuts[1] && pa[2]>pos[2]-cuts[2] && pa[2]<pos[2]+cuts[2]){ | |
198 | distance = (TMath::Sqrt(pow((pa[0]-pb[0]),2) + pow((pa[1]-pb[1]),2) + pow((pa[2]-pb[2]),2))); | |
199 | if(distance<distancecut){ | |
200 | hx->Fill(pa[0]); hy->Fill(pa[1]); hz->Fill(pa[2]); | |
201 | hx->Fill(pb[0]); hy->Fill(pb[1]); hz->Fill(pb[2]); | |
202 | ncontr++; | |
203 | } | |
204 | } | |
205 | } | |
206 | ||
207 | // Third line is formed by p1-p4 and fourth line by p3-p2 | |
208 | // We find two points on each line which form the closest distance of the two lines | |
209 | // pa[0],pa[1],pa[2]: points on line 3 and pb[0],pb[1],pb[2]: points on line 4 | |
210 | // Next: Consider x, y and z to be less than cuts[0], cuts[1] and cuts[2], respectively | |
211 | if(TMath::Abs(hphi1-hphi4)<phicut && TMath::Abs(hphi3-hphi2)<phicut) { | |
212 | CalculateLine(p1, p4, p3, p2, pa, pb); | |
213 | if (pa[0]>pos[0]-cuts[0] && pa[0]<pos[0]+cuts[0] && pa[1]>pos[1]-cuts[1] && pa[1]<pos[1]+cuts[1]){ | |
214 | distance = (TMath::Sqrt(pow((pa[0]-pb[0]),2) + pow((pa[1]-pb[1]),2) + pow((pa[2]-pb[2]),2))); | |
215 | if(distance<distancecut){ | |
216 | hx->Fill(pa[0]); hy->Fill(pa[1]); hz->Fill(pa[2]); | |
217 | hx->Fill(pb[0]); hy->Fill(pb[1]); hz->Fill(pb[2]); | |
218 | ncontr++; | |
219 | } | |
220 | } | |
221 | } | |
222 | } | |
223 | } | |
224 | ||
225 | Int_t maxbinx = hx->GetMaximumBin(); | |
226 | Int_t maxbiny = hy->GetMaximumBin(); | |
227 | Int_t maxbinz = hz->GetMaximumBin(); | |
228 | pos[0] = ilow + ((ihigh-ilow)/ibin)*maxbinx; | |
229 | pos[1] = ilow + ((ihigh-ilow)/ibin)*maxbiny; | |
230 | pos[2] = ilowz + ((ihighz-ilowz)/ibinz)*maxbinz; | |
231 | hx->Delete(); | |
232 | hy->Delete(); | |
233 | hz->Delete(); | |
234 | } | |
235 | ||
236 | void AliITSVertexer3DTapan::CalculateVertex3d2(Double_t pos[3], Float_t cuts[3], Int_t &ncontr, Double_t sigpos[3]) { | |
237 | // | |
238 | // This function reconstructs second two steps of vertex | |
239 | // | |
240 | ||
241 | Double_t p1[4], p2[4], p3[4], p4[4]; | |
242 | Double_t pa[3], pb[3]; | |
243 | Double_t hphi1, hphi2, hphi3, hphi4; | |
244 | ||
245 | ncontr = 0; | |
246 | Float_t phicut = 0.3; | |
247 | Double_t distance; Float_t distancecut = 1.0; | |
248 | ||
249 | Double_t vertx =0.; Double_t verty =0.; Double_t vertz =0.; | |
250 | Double_t vertx2 =0.; Double_t verty2 =0.; Double_t vertz2 =0.; | |
251 | ||
252 | for (Int_t ip1=0; ip1<ficlu1; ip1++) { | |
253 | // Two points on layer1: p1 and p3 | |
254 | p1[0] = fX1[ip1]; p1[1] = fY1[ip1]; p1[2] = fZ1[ip1]; | |
255 | p3[0] = fX1[ip1+1]; p3[1] = fY1[ip1+1]; p3[2] = fZ1[ip1+1]; | |
256 | hphi1 = fPhi1[ip1]; hphi3 = fPhi1[ip1+1]; | |
257 | ||
258 | for (Int_t ip2=0; ip2<ficlu2; ip2++) { | |
259 | // Two points on layer 2: p2 and p4 | |
260 | p2[0] = fX2[ip2]; p2[1] = fY2[ip2]; p2[2] = fZ2[ip2]; | |
261 | p4[0] = fX2[ip2+1]; p4[1] = fY2[ip2+1]; p4[2] = fZ2[ip2+1]; | |
262 | hphi2 = fPhi2[ip2]; hphi4 = fPhi2[ip2+1]; | |
263 | ||
264 | // First line is formed by p1-p2 and second line by p3-p4 | |
265 | // We find two points on each line which form the closest distance of the two lines | |
266 | // pa[0],pa[1],pa[2] are the points on line 1 and pb[0],pb[1],pb[2] are the points on line 2 | |
267 | // Next: Consider x, y and z to be less than cuts[0], cuts[1] and cuts[2], respectively | |
268 | ||
269 | if(TMath::Abs(hphi1-hphi2)<phicut && TMath::Abs(hphi3-hphi4)<phicut){ | |
270 | CalculateLine(p1, p2, p3, p4, pa, pb); | |
271 | ||
272 | // We consider the points where x, y and z points are less than xcut, ycut and zcut, respectively | |
273 | if (pa[0]>pos[0]-cuts[0] && pa[0]<pos[0]+cuts[0] && pa[1]>pos[1]-cuts[1] && pa[1]<pos[1]+cuts[1] && pa[2]>pos[2]-cuts[2] && pa[2]<pos[2]+cuts[2]){ | |
274 | distance = (TMath::Sqrt(pow((pa[0]-pb[0]),2) + pow((pa[1]-pb[1]),2) + pow((pa[2]-pb[2]),2))); | |
275 | if(distance<distancecut){ | |
276 | ncontr++; | |
277 | vertx = vertx + pa[0]; verty = verty + pa[1]; vertz = vertz + pa[2]; | |
278 | vertx2 = vertx2 + pa[0]*pa[0]; verty2 = verty2 + pa[1]*pa[1]; vertz2 = vertz2 + pa[2]*pa[2]; | |
279 | ncontr++; | |
280 | vertx = vertx + pb[0]; verty = verty + pb[1]; vertz = vertz + pb[2]; | |
281 | vertx2 = vertx2 + pb[0]*pb[0]; verty2 = verty2 + pb[1]*pb[1]; vertz2 = vertz2 + pb[2]*pb[2]; | |
282 | } | |
283 | } | |
284 | } | |
285 | ||
286 | // Third line is formed by p1-p4 and fourth line by p3-p2 | |
287 | // We find two points on each line which form the closest distance of the two lines | |
288 | // pa[0],pa[1],pa[2] are the points on line 3 and pb[0],pb[1],pb[2] are the points on line 4 | |
289 | // Next: Consider x, y and z to be less than cuts[0], cuts[1] and cuts[2], respectively | |
290 | if(TMath::Abs(hphi1-hphi4)<phicut && TMath::Abs(hphi3-hphi2)<phicut) { | |
291 | ||
292 | CalculateLine(p1, p4, p3, p2, pa, pb); | |
293 | if (pa[0]>pos[0]-cuts[0] && pa[0]<pos[0]+cuts[0] && pa[1]>pos[1]-cuts[1] && pa[1]<pos[1]+cuts[1] && pa[2]>pos[2]-cuts[2] && pa[2]<pos[2]+cuts[2]){ | |
294 | distance = (TMath::Sqrt(pow((pa[0]-pb[0]),2) + pow((pa[1]-pb[1]),2) + pow((pa[2]-pb[2]),2))); | |
295 | if(distance<distancecut){ | |
296 | ncontr++; | |
297 | vertx = vertx + pa[0]; verty = verty + pa[1]; vertz = vertz + pa[2]; | |
298 | vertx2 = vertx2 + pa[0]*pa[0]; verty2 = verty2 + pa[1]*pa[1]; vertz2 = vertz2 + pa[2]*pa[2]; | |
299 | ncontr++; | |
300 | vertx = vertx + pb[0]; verty = verty + pb[1]; vertz = vertz + pb[2]; | |
301 | vertx2 = vertx2 + pb[0]*pb[0]; verty2 = verty2 + pb[1]*pb[1]; vertz2 = vertz2 + pb[2]*pb[2]; | |
302 | } | |
303 | } | |
304 | } | |
305 | } | |
306 | } | |
307 | ||
308 | if(ncontr>0){ | |
309 | pos[0] = vertx/ncontr; pos[1] = verty/ncontr; pos[2] = vertz/ncontr; | |
310 | vertx2 = vertx2/ncontr; verty2 = verty2/ncontr; vertz2 = vertz2/ncontr; | |
311 | sigpos[0] = TMath::Sqrt(vertx2 - pos[0]*pos[0]); | |
312 | sigpos[1] = TMath::Sqrt(verty2 - pos[1]*pos[1]); | |
313 | sigpos[2] = TMath::Sqrt(vertz2 - pos[2]*pos[2]); | |
314 | } | |
315 | ncontr = ncontr/2; | |
316 | } | |
317 | ||
318 | void AliITSVertexer3DTapan::CalculatePhi(Float_t fx, Float_t fy, Float_t & phi) | |
319 | { | |
320 | //calculates phi | |
321 | Float_t ybyx, phi1; | |
322 | const Float_t kradian = 180./3.141592654; | |
323 | ||
324 | if(fx==0.) | |
325 | { | |
326 | if(fy>0.) phi = 90.; | |
327 | if(fy<0.) phi = 270.; | |
328 | } | |
329 | if(fx != 0.) | |
330 | { | |
331 | ybyx = fy/fx; | |
332 | if(ybyx < 0) ybyx = - ybyx; | |
333 | phi1 = TMath::ATan(ybyx)*kradian; | |
334 | if(fx > 0 && fy > 0) phi = phi1; // 1st Quadrant | |
335 | if(fx < 0 && fy > 0) phi = 180 - phi1; // 2nd Quadrant | |
336 | if(fx < 0 && fy < 0) phi = 180 + phi1; // 3rd Quadrant | |
337 | if(fx > 0 && fy < 0) phi = 360 - phi1; // 4th Quadrant | |
338 | ||
339 | } | |
340 | phi = phi/kradian; | |
341 | } | |
342 | ||
343 | void AliITSVertexer3DTapan::CalculateLine(Double_t p1[4], Double_t p2[4], Double_t p3[4], Double_t p4[4], Double_t pa[3], Double_t pb[3]) const{ | |
344 | //calculates line | |
345 | Double_t p13x, p13y, p13z; | |
346 | Double_t p21x, p21y, p21z; | |
347 | Double_t p43x, p43y, p43z; | |
348 | Double_t d1343, d4321, d1321, d4343, d2121; | |
349 | Double_t numer, denom; | |
350 | Double_t mua, mub; | |
351 | mua = 0.; mub = 0.; | |
352 | ||
353 | p13x = p1[0] - p3[0]; | |
354 | p13y = p1[1] - p3[1]; | |
355 | p13z = p1[2] - p3[2]; | |
356 | ||
357 | p21x = p2[0] - p1[0]; | |
358 | p21y = p2[1] - p1[1]; | |
359 | p21z = p2[2] - p1[2]; | |
360 | ||
361 | p43x = p4[0] - p3[0]; | |
362 | p43y = p4[1] - p3[1]; | |
363 | p43z = p4[2] - p3[2]; | |
364 | ||
365 | d1343 = p13x * p43x + p13y * p43y + p13z * p43z; | |
366 | d4321 = p43x * p21x + p43y * p21y + p43z * p21z; | |
367 | d1321 = p13x * p21x + p13y * p21y + p13z * p21z; | |
368 | d4343 = p43x * p43x + p43y * p43y + p43z * p43z; | |
369 | d2121 = p21x * p21x + p21y * p21y + p21z * p21z; | |
370 | ||
371 | denom = d2121 * d4343 - d4321 * d4321; | |
372 | numer = d1343 * d4321 - d1321 * d4343; | |
373 | ||
374 | if(denom>0) mua = numer / denom; | |
375 | if(d4343>0) mub = (d1343 + d4321 * (mua)) / d4343; | |
376 | ||
377 | pa[0] = p1[0] + mua * p21x; | |
378 | pa[1] = p1[1] + mua * p21y; | |
379 | pa[2] = p1[2] + mua * p21z; | |
380 | ||
381 | pb[0] = p3[0] + mub * p43x; | |
382 | pb[1] = p3[1] + mub * p43y; | |
383 | pb[2] = p3[2] + mub * p43z; | |
384 | } | |
385 |