Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, 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 : // Implementation of the cascade vertexer class
18 : // Reads V0s and tracks, writes out cascade vertices
19 : // Fills the ESD with the cascades
20 : // Origin: Christian Kuhn, IReS, Strasbourg, christian.kuhn@ires.in2p3.fr
21 : //-------------------------------------------------------------------------
22 :
23 : //modified by R. Vernet 30/6/2006 : daughter label
24 : //modified by R. Vernet 3/7/2006 : causality
25 : //modified by I. Belikov 24/11/2006 : static setter for the default cuts
26 :
27 : #include "AliESDEvent.h"
28 : #include "AliESDcascade.h"
29 : #include "AliCascadeVertexer.h"
30 :
31 172 : ClassImp(AliCascadeVertexer)
32 :
33 : //A set of loose cuts
34 : Double_t
35 : AliCascadeVertexer::fgChi2max=33.; //maximal allowed chi2
36 : Double_t
37 : AliCascadeVertexer::fgDV0min=0.01; //min V0 impact parameter
38 : Double_t
39 : AliCascadeVertexer::fgMassWin=0.008; //"window" around the Lambda mass
40 : Double_t
41 : AliCascadeVertexer::fgDBachMin=0.01; //min bachelor impact parameter
42 : Double_t
43 : AliCascadeVertexer::fgDCAmax=2.0; //max DCA between the V0 and the track
44 : Double_t
45 : AliCascadeVertexer::fgCPAmin=0.98; //min cosine of the cascade pointing angle
46 : Double_t
47 : AliCascadeVertexer::fgRmin=0.2; //min radius of the fiducial volume
48 : Double_t
49 : AliCascadeVertexer::fgRmax=100.; //max radius of the fiducial volume
50 :
51 :
52 : Int_t AliCascadeVertexer::V0sTracks2CascadeVertices(AliESDEvent *event) {
53 : //--------------------------------------------------------------------
54 : // This function reconstructs cascade vertices
55 : // Adapted to the ESD by I.Belikov (Jouri.Belikov@cern.ch)
56 : //--------------------------------------------------------------------
57 16 : const AliESDVertex *vtxT3D=event->GetPrimaryVertex();
58 :
59 8 : Double_t xPrimaryVertex=vtxT3D->GetX();
60 8 : Double_t yPrimaryVertex=vtxT3D->GetY();
61 8 : Double_t zPrimaryVertex=vtxT3D->GetZ();
62 :
63 8 : Double_t b=event->GetMagneticField();
64 8 : Int_t nV0=(Int_t)event->GetNumberOfV0s();
65 :
66 : //stores relevant V0s in an array
67 8 : TObjArray vtcs(nV0);
68 : Int_t i;
69 72 : for (i=0; i<nV0; i++) {
70 28 : AliESDv0 *v=event->GetV0(i);
71 42 : if (v->GetOnFlyStatus()) continue;
72 28 : if (v->GetD(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex)<fDV0min) continue;
73 14 : vtcs.AddLast(v);
74 14 : }
75 8 : nV0=vtcs.GetEntriesFast();
76 :
77 : // stores relevant tracks in another array
78 8 : Int_t nentr=(Int_t)event->GetNumberOfTracks();
79 8 : int trk[nentr], ntr=0;
80 320 : for (i=0; i<nentr; i++) {
81 152 : AliESDtrack *esdtr=event->GetTrack(i);
82 152 : ULong_t status=esdtr->GetStatus();
83 152 : if (status&AliESDtrack::kITSpureSA) continue;
84 152 : if ((status&AliESDtrack::kITSrefit)==0)
85 60 : if ((status&AliESDtrack::kTPCrefit)==0) continue;
86 :
87 376 : if (TMath::Abs(esdtr->GetD(xPrimaryVertex,yPrimaryVertex,b))<fDBachMin) continue;
88 :
89 80 : trk[ntr++]=i;
90 80 : }
91 :
92 : Double_t massLambda=1.11568;
93 : Int_t ncasc=0;
94 :
95 : // Looking for the cascades...
96 :
97 44 : for (i=0; i<nV0; i++) { //loop on V0s
98 :
99 14 : AliESDv0 *v=(AliESDv0*)vtcs.UncheckedAt(i);
100 14 : AliESDv0 v0(*v);
101 14 : v0.ChangeMassHypothesis(kLambda0); // the v0 must be Lambda
102 26 : if (TMath::Abs(v0.GetEffMass()-massLambda)>fMassWin) continue;
103 :
104 28 : for (Int_t j=0; j<ntr; j++) {//loop on tracks
105 12 : Int_t bidx=trk[j];
106 : //Bo: if (bidx==v->GetNindex()) continue; //bachelor and v0's negative tracks must be different
107 14 : if (bidx==v0.GetIndex(0)) continue; //Bo: consistency 0 for neg
108 10 : AliESDtrack *btrk=event->GetTrack(bidx);
109 26 : if (btrk->GetSign()>0) continue; // bachelor's charge
110 :
111 : AliESDv0 *pv0=&v0;
112 4 : AliExternalTrackParam bt(*btrk), *pbt=&bt;
113 :
114 4 : Double_t dca=PropagateToDCA(pv0,pbt,b);
115 4 : if (dca > fDCAmax) continue;
116 :
117 4 : AliESDcascade cascade(*pv0,*pbt,bidx);//constucts a cascade candidate
118 : //PH if (cascade.GetChi2Xi() > fChi2max) continue;
119 :
120 4 : Double_t x,y,z; cascade.GetXYZcascade(x,y,z); // Bo: bug correction
121 4 : Double_t r2=x*x + y*y;
122 6 : if (r2 > fRmax2) continue; // condition on fiducial zone
123 2 : if (r2 < fRmin2) continue;
124 :
125 2 : Double_t pxV0,pyV0,pzV0;
126 2 : pv0->GetPxPyPz(pxV0,pyV0,pzV0);
127 2 : if (x*pxV0+y*pyV0+z*pzV0 < 0) continue; //causality
128 :
129 2 : Double_t x1,y1,z1; pv0->GetXYZ(x1,y1,z1);
130 2 : if (r2 > (x1*x1+y1*y1)) continue;
131 :
132 6 : if (cascade.GetCascadeCosineOfPointingAngle(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex) <fCPAmin) continue; //condition on the cascade pointing angle
133 :
134 0 : cascade.SetDcaXiDaughters(dca);
135 0 : event->AddCascade(&cascade);
136 0 : ncasc++;
137 12 : } // end loop tracks
138 16 : } // end loop V0s
139 :
140 : // Looking for the anti-cascades...
141 :
142 44 : for (i=0; i<nV0; i++) { //loop on V0s
143 14 : AliESDv0 *v=(AliESDv0*)vtcs.UncheckedAt(i);
144 14 : AliESDv0 v0(*v);
145 14 : v0.ChangeMassHypothesis(kLambda0Bar); //the v0 must be anti-Lambda
146 28 : if (TMath::Abs(v0.GetEffMass()-massLambda)>fMassWin) continue;
147 :
148 0 : for (Int_t j=0; j<ntr; j++) {//loop on tracks
149 0 : Int_t bidx=trk[j];
150 : //Bo: if (bidx==v->GetPindex()) continue; //bachelor and v0's positive tracks must be different
151 0 : if (bidx==v0.GetIndex(1)) continue; //Bo: consistency 1 for pos
152 0 : AliESDtrack *btrk=event->GetTrack(bidx);
153 0 : if (btrk->GetSign()<0) continue; // bachelor's charge
154 :
155 : AliESDv0 *pv0=&v0;
156 0 : AliExternalTrackParam bt(*btrk), *pbt=&bt;
157 :
158 0 : Double_t dca=PropagateToDCA(pv0,pbt,b);
159 0 : if (dca > fDCAmax) continue;
160 :
161 0 : AliESDcascade cascade(*pv0,*pbt,bidx); //constucts a cascade candidate
162 : //PH if (cascade.GetChi2Xi() > fChi2max) continue;
163 :
164 0 : Double_t x,y,z; cascade.GetXYZcascade(x,y,z); // Bo: bug correction
165 0 : Double_t r2=x*x + y*y;
166 0 : if (r2 > fRmax2) continue; // condition on fiducial zone
167 0 : if (r2 < fRmin2) continue;
168 :
169 0 : Double_t pxV0,pyV0,pzV0;
170 0 : pv0->GetPxPyPz(pxV0,pyV0,pzV0);
171 0 : if (x*pxV0+y*pyV0+z*pzV0 < 0) continue; //causality
172 :
173 0 : Double_t x1,y1,z1; pv0->GetXYZ(x1,y1,z1);
174 0 : if (r2 > (x1*x1+y1*y1)) continue;
175 :
176 0 : if (cascade.GetCascadeCosineOfPointingAngle(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex) < fCPAmin) continue; //condition on the cascade pointing angle
177 :
178 0 : cascade.SetDcaXiDaughters(dca);
179 0 : event->AddCascade(&cascade);
180 0 : ncasc++;
181 :
182 0 : } // end loop tracks
183 14 : } // end loop V0s
184 :
185 8 : Info("V0sTracks2CascadeVertices","Number of reconstructed cascades: %d",ncasc);
186 :
187 : return 0;
188 8 : }
189 :
190 :
191 : Double_t AliCascadeVertexer::Det(Double_t a00, Double_t a01, Double_t a10, Double_t a11) const {
192 : //--------------------------------------------------------------------
193 : // This function calculates locally a 2x2 determinant
194 : //--------------------------------------------------------------------
195 96 : return a00*a11 - a01*a10;
196 : }
197 :
198 : Double_t AliCascadeVertexer::Det(Double_t a00,Double_t a01,Double_t a02,
199 : Double_t a10,Double_t a11,Double_t a12,
200 : Double_t a20,Double_t a21,Double_t a22) const {
201 : //--------------------------------------------------------------------
202 : // This function calculates locally a 3x3 determinant
203 : //--------------------------------------------------------------------
204 24 : return a00*Det(a11,a12,a21,a22)-a01*Det(a10,a12,a20,a22)+a02*Det(a10,a11,a20,a21);
205 : }
206 :
207 :
208 :
209 :
210 : Double_t AliCascadeVertexer::PropagateToDCA(AliESDv0 *v, AliExternalTrackParam *t, Double_t b) {
211 : //--------------------------------------------------------------------
212 : // This function returns the DCA between the V0 and the track
213 : //--------------------------------------------------------------------
214 8 : Double_t alpha=t->GetAlpha(), cs1=TMath::Cos(alpha), sn1=TMath::Sin(alpha);
215 4 : Double_t r[3]; t->GetXYZ(r);
216 4 : Double_t x1=r[0], y1=r[1], z1=r[2];
217 4 : Double_t p[3]; t->GetPxPyPz(p);
218 4 : Double_t px1=p[0], py1=p[1], pz1=p[2];
219 :
220 4 : Double_t x2,y2,z2; // position and momentum of V0
221 4 : Double_t px2,py2,pz2;
222 :
223 4 : v->GetXYZ(x2,y2,z2);
224 4 : v->GetPxPyPz(px2,py2,pz2);
225 :
226 : // calculation dca
227 :
228 4 : Double_t dd= Det(x2-x1,y2-y1,z2-z1,px1,py1,pz1,px2,py2,pz2);
229 4 : Double_t ax= Det(py1,pz1,py2,pz2);
230 4 : Double_t ay=-Det(px1,pz1,px2,pz2);
231 4 : Double_t az= Det(px1,py1,px2,py2);
232 :
233 4 : Double_t dca=TMath::Abs(dd)/TMath::Sqrt(ax*ax + ay*ay + az*az);
234 4 : if (dca > fDCAmax) return 1.e+33;
235 :
236 : //points of the DCA
237 8 : Double_t t1 = Det(x2-x1,y2-y1,z2-z1,px2,py2,pz2,ax,ay,az)/
238 4 : Det(px1,py1,pz1,px2,py2,pz2,ax,ay,az);
239 :
240 4 : x1 += px1*t1; y1 += py1*t1; //z1 += pz1*t1;
241 :
242 4 : if (x1*x1+y1*y1 > fRmaxMargin2) return 1.e+33;
243 :
244 : //propagate track to the points of DCA
245 :
246 4 : x1=x1*cs1 + y1*sn1;
247 :
248 4 : if (!t->PropagateTo(x1,b)) {
249 0 : AliError("Propagation failed");
250 : // AliErrorF("Propagation failed for X=%f | V0: %f %f %f",x1,x2,y2,z2);
251 : // t->Print();
252 : //
253 0 : return 1.e+33;
254 : }
255 :
256 4 : return dca;
257 4 : }
258 :
259 :
260 :
261 :
262 :
263 :
264 :
265 :
266 :
267 :
268 :
|
