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