]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSVertexer3DTapan.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
HLT mode retrieved from the OCDB and used to set SDD raw data format in simulation...
[u/mrichter/AliRoot.git] / ITS / AliITSVertexer3DTapan.cxx
CommitLineData
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
31ClassImp(AliITSVertexer3DTapan)
32
308c2f7c 33void 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 97AliESDVertex *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
158void 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
236void 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
318void 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
343void 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