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