Added missing files and defines, compilation should now work again.
[u/mrichter/AliRoot.git] / HLT / trigger / AliD0Trigger.cxx
CommitLineData
5a31e9df 1#include "AliD0Trigger.h"
2#include <TMath.h>
3#include "AliITStrackV2.h"
4#include "AliL3Transform.h"
5#include <TVector3.h>
0bd0c1ef 6#include <iostream.h>
5a31e9df 7
8ClassImp(AliD0Trigger)
9
10AliD0Trigger::AliD0Trigger()
11{
12 posTrack=0;
13 negTrack=0;
14}
0bd0c1ef 15AliD0Trigger::AliD0Trigger(double c[6],double b,double pv[3])
5a31e9df 16{
17 posTrack=0;
18 negTrack=0;
19 cutV0low=c[0];
20 cutV0high=c[1];
21 cutInvMass=c[2];
22 cutPointAngle=c[3];
0bd0c1ef 23 cutd0d0=c[4];
24 cutCosThetaStar=c[5];
25 cutpTchild=c[6];
5a31e9df 26 Bfield=b;
27 primaryVertex[0]=pv[0];
28 primaryVertex[1]=pv[1];
29 primaryVertex[2]=pv[2];
30}
31AliD0Trigger::AliD0Trigger(AliITStrackV2 * posT, AliITStrackV2 * negT)
32{
33 posTrack=posT;
34 negTrack=negT;
35}
36
37
38AliD0Trigger::~AliD0Trigger()
39{
40
41}
42
43void AliD0Trigger::SetTracks(AliITStrackV2 * posT, AliITStrackV2 * negT){
44 posTrack=posT;
45 negTrack=negT;
46}
47bool AliD0Trigger::FindInvMass() {
48
49 double mD0 = 1.8645;
50 double mPi = 0.13957;
51 double mKa = 0.49368;
52
53 double mom2[2],momTot2,energy[2];
54
55 mom2[0]=pow(momenta[0],2)+pow(momenta[1],2)+pow(momenta[2],2);
56 mom2[1]=pow(momenta[3],2)+pow(momenta[4],2)+pow(momenta[5],2);
57
58 momTot2 = pow(momenta[0]+momenta[3],2)+pow(momenta[1]+momenta[4],2)+pow(momenta[2]+momenta[5],2);
59
60 // D0 -> K- Pi+
61 energy[1] = sqrt(pow(mKa,2)+mom2[1]);
62 energy[0] = sqrt(pow(mPi,2)+mom2[0]);
63
64 double minvD0 = sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
65
66 // D0bar -> K+ Pi-
67 energy[0] = sqrt(mKa*mKa+mom2[0]);
68 energy[1] = sqrt(mPi*mPi+mom2[1]);
69
70 double minvD0bar = sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
71
72 if(fabs(minvD0-mD0) < cutInvMass) return true;
73 if(fabs(minvD0bar-mD0) < cutInvMass) return true;
74 return false;
75
76}
77bool AliD0Trigger::FindV0(){
78
0bd0c1ef 79 bool goodV0=false;
80
81 double Gxpos=0,Gypos=0,Gxneg=0,Gyneg=0;
82 double r1=0,r2=0,a1=0,a2=0,b1=0,b2=0;
83 double q=0,p=0,t=0,Fa=0,Fb=0,Fc=0;
84 double y1=0,y2=0,x1=0,x2=0,x3=0,x4=0;
85 double V01=0,V02=0,V03=0,V04=0;
86 double V0[4]={0,0,0,0};
87 bestV0[0]=0; bestV0[1]=0; bestV0[2]=0;
5a31e9df 88
0bd0c1ef 89 Gxpos=posTrack->GetX()*cos(posTrack->GetAlpha())-posTrack->GetY()*sin(posTrack->GetAlpha());
90 Gypos=posTrack->GetX()*sin(posTrack->GetAlpha())+posTrack->GetY()*cos(posTrack->GetAlpha());
91 Gxneg=negTrack->GetX()*cos(negTrack->GetAlpha())-negTrack->GetY()*sin(negTrack->GetAlpha());
92 Gyneg=negTrack->GetX()*sin(negTrack->GetAlpha())+negTrack->GetY()*cos(negTrack->GetAlpha());
5a31e9df 93
0bd0c1ef 94 r1=fabs(1/(AliL3Transform::GetBFact()*Bfield*posTrack->Get1Pt()));
95 r2=fabs(1/(AliL3Transform::GetBFact()*Bfield*negTrack->Get1Pt()));
96
97 a1=Gxpos-(r1*cos(posTrack->GetAlpha()+asin(posTrack->GetSnp())+TMath::PiOver2()));
98 a2=Gxneg-(r2*cos(negTrack->GetAlpha()+asin(negTrack->GetSnp())-TMath::PiOver2()));
99 b1=Gypos-(r1*sin(posTrack->GetAlpha()+asin(posTrack->GetSnp())+TMath::PiOver2()));
100 b2=Gyneg-(r2*sin(negTrack->GetAlpha()+asin(negTrack->GetSnp())-TMath::PiOver2()));
101
102 //double a3=a1-a2;
103 //double b3=b1-b2;
104 //double r3=sqrt(pow(a3,2)+pow(b3,2));
105
106 //if(r3<r1+r2){
107 q=a1-a2;
108 p=2*(b1-b2);
109 t=pow(r2,2)-pow(r2,2)+pow(b1,2)-pow(b2,2)-pow(q,2);
110
111 Fa=pow(p,2)+(4*pow(q,2));
112 Fb=-(2*t*p+8*pow(q,2)*b1);
113 Fc=pow(t,2)-(4*pow(q,2)*pow(r1,2))+(4*pow(q,2)*pow(b1,2));
114
115 if(pow(Fb,2)-(4*Fa*Fc)>=0){
116 y1=(-Fb+(sqrt(pow(Fb,2)-(4*Fa*Fc))))/(2*Fa); //noe feil her. floating point
117 y2=(-Fb-(sqrt(pow(Fb,2)-(4*Fa*Fc))))/(2*Fa); //trolig negativ under rot
118
119 x1=sqrt(pow(r1,2)-pow((y1-b1),2))+a1;
120 x2=sqrt(pow(r1,2)-pow((y2-b1),2))+a1;
121 x3=-sqrt(pow(r1,2)-pow((y1-b1),2))+a1;
122 x4=-sqrt(pow(r1,2)-pow((y2-b1),2))+a1;
123
124 V01=sqrt(pow(x1-primaryVertex[0],2)+pow(y1-primaryVertex[1],2));
125 V02=sqrt(pow(x2-primaryVertex[0],2)+pow(y2-primaryVertex[1],2));
126 V03=sqrt(pow(x3-primaryVertex[0],2)+pow(y1-primaryVertex[1],2));
127 V04=sqrt(pow(x4-primaryVertex[0],2)+pow(y2-primaryVertex[1],2));
128
129 V0[0]=V01;V0[1]=V02;V0[2]=V03;V0[3]=V04;
130
131 int index=0;
132
133 double nearestV0=V0[0];
134
135 for(int i=0;i<4;i++){
136 if(nearestV0>V0[i]){
137 nearestV0=V0[i];
138 index=i;
5a31e9df 139 }
0bd0c1ef 140 }
141 cout<<"index: "<<index<<endl;
142 if(nearestV0<cutV0high){
143 if(nearestV0>cutV0low){
144 goodV0=true;
145 if(index==0){
146 bestV0[0]=x1; bestV0[1]=y1; bestV0[2]=0.;
147 }
148 if(index==1){
149 bestV0[0]=x2; bestV0[1]=y2; bestV0[2]=0.;
150 }
151 if(index==2){
152 bestV0[0]=x3; bestV0[1]=y1; bestV0[2]=0.;
153 }
154 if(index==3){
155 bestV0[0]=x4; bestV0[1]=y2; bestV0[2]=0.;
5a31e9df 156 }
157 }
158 }
0bd0c1ef 159 //}
5a31e9df 160 }
0bd0c1ef 161
162 cout<<"My V0: x: "<<bestV0[0]<<" y: "<<bestV0[1]<<" z: "<<bestV0[2]<<endl;
163
5a31e9df 164 return goodV0;
165}
166void AliD0Trigger::FindMomentaAtVertex(){
167
168 double r1=fabs(1/(AliL3Transform::GetBFact()*Bfield*posTrack->Get1Pt()));
169 double r2=fabs(1/(AliL3Transform::GetBFact()*Bfield*negTrack->Get1Pt()));
170
171 double Gx1=posTrack->GetX()*cos(posTrack->GetAlpha())-posTrack->GetY()*sin(posTrack->GetAlpha());
172 double Gy1=posTrack->GetX()*sin(posTrack->GetAlpha())+posTrack->GetY()*cos(posTrack->GetAlpha());
173 double Gx2=negTrack->GetX()*cos(negTrack->GetAlpha())-negTrack->GetY()*sin(negTrack->GetAlpha());
174 double Gy2=negTrack->GetX()*sin(negTrack->GetAlpha())+negTrack->GetY()*cos(negTrack->GetAlpha());
175
176 double centerx1=Gx1-(r1*cos(posTrack->GetAlpha()+asin(posTrack->GetSnp())+TMath::PiOver2()));
177 double centerx2=Gx2-(r2*cos(negTrack->GetAlpha()+asin(negTrack->GetSnp())-TMath::PiOver2()));
178 double centery1=Gy1-(r1*sin(posTrack->GetAlpha()+asin(posTrack->GetSnp())+TMath::PiOver2()));
179 double centery2=Gy2-(r2*sin(negTrack->GetAlpha()+asin(negTrack->GetSnp())-TMath::PiOver2()));
180
181 double a1=sqrt(pow(Gx1-bestV0[0],2)+pow(Gy1-bestV0[1],2));
182 double b1=sqrt(pow(bestV0[0]-centerx1,2)+pow(bestV0[1]-centery1,2));
183 double c1=sqrt(pow(Gx1-centerx1,2)+pow(Gy1-centery1,2));
184 double a2=sqrt(pow(Gx2-bestV0[0],2)+pow(Gy2-bestV0[1],2));
185 double b2=sqrt(pow(bestV0[0]-centerx2,2)+pow(bestV0[1]-centery2,2));
186 double c2=sqrt(pow(Gx2-centerx2,2)+pow(Gy2-centery2,2));
187
188 double alpha1=acos((pow(b1,2)+pow(c1,2)-pow(a1,2))/(2*b1*c1));
189 double alpha2=acos((pow(b2,2)+pow(c2,2)-pow(a2,2))/(2*b2*c2));
190
191 double sign1=0;
192 double sign2=0;
193
194 if((sqrt(pow(Gx1,2)+pow(Gy1,2)))>(sqrt(pow(bestV0[0],2)+pow(bestV0[1],2))))
195 sign1=-1;
196 else
197 sign1=1;
198 if((sqrt(pow(Gx2,2)+pow(Gy2,2)))>(sqrt(pow(bestV0[0],2)+pow(bestV0[1],2))))
199 sign2=1;
200 else
201 sign2=-1;
202
203 double psi1=(posTrack->GetAlpha()+asin(posTrack->GetSnp()))+(sign1*alpha1);
204 double psi2=(negTrack->GetAlpha()+asin(negTrack->GetSnp()))+(sign2*alpha2);
205
206 double ptP = 1./fabs(posTrack->Get1Pt());
207 momenta[0] = ptP*cos(psi1);
208 momenta[1] = ptP*sin(psi1);
209 momenta[2] = ptP*posTrack->GetTgl();
210
211 double ptN = 1./fabs(negTrack->Get1Pt());
212 momenta[3] = ptN*cos(psi2);
213 momenta[4] = ptN*sin(psi2);
214 momenta[5] = ptN*negTrack->GetTgl();
215
216}
217Bool_t AliD0Trigger::PointingAngle(){
218
219 TVector3 mom(momenta[0]+momenta[3],momenta[1]+momenta[4],momenta[2]+momenta[5]);
220 TVector3 flight(bestV0[0]-primaryVertex[0],bestV0[1]-primaryVertex[1],bestV0[2]-primaryVertex[2]);
221
222 double pta = mom.Angle(flight);
223
224 if(cos(pta)<cutPointAngle)
225 return true;
226 else
227 return false;
228}
0bd0c1ef 229void AliD0Trigger::SetMomenta(double m[6])
230{
231 momenta[0] = m[0];
232 momenta[1] = m[1];
233 momenta[2] = m[2];
234 momenta[3] = m[3];
235 momenta[4] = m[4];
236 momenta[5] = m[5];
237}
238void AliD0Trigger::FindMomentaOffline()
239{
240 double ptP,alphaP,phiP,ptN,alphaN,phiN;
241 // momenta of the tracks at the vertex
242 ptP = 1./TMath::Abs(posTrack->Get1Pt());
243 alphaP = posTrack->GetAlpha();
244 phiP = alphaP+TMath::ASin(posTrack->GetSnp());
245 momenta[0] = ptP*TMath::Cos(phiP);
246 momenta[1] = ptP*TMath::Sin(phiP);
247 momenta[2] = ptP*posTrack->GetTgl();
248
249 ptN = 1./TMath::Abs(negTrack->Get1Pt());
250 alphaN = negTrack->GetAlpha();
251 phiN = alphaN+TMath::ASin(negTrack->GetSnp());
252 momenta[3] = ptN*TMath::Cos(phiN);
253 momenta[4] = ptN*TMath::Sin(phiN);
254 momenta[5] = ptN*negTrack->GetTgl();
255}
256bool AliD0Trigger::FindV0offline(double v[3])
257{
258
259 bool goodV0=false;
260 double r=0;
261
262 bestV0[0]=v[0];
263 bestV0[1]=v[1];
264 bestV0[2]=v[2];
265
266 r=sqrt(pow(bestV0[0]-primaryVertex[0],2)+pow(bestV0[1]-primaryVertex[1],2)+pow(bestV0[2]-primaryVertex[2],2));
267
268 if(r<cutV0high){
269 if(r>cutV0low){
270 goodV0=true;
271 }
272 }
273
274 return goodV0;
275
276}
277bool AliD0Trigger::d0d0()
278{
279 bool goodd0=false;
280 double d00=0, d01=0;
281
282 d00 = 10000.*posTrack->GetD(primaryVertex[0],primaryVertex[1]);
283 d01 = -10000.*negTrack->GetD(primaryVertex[0],primaryVertex[1]);
284
285 if(d00*d01<cutd0d0){
286 goodd0=true;
287 }
288
289 return goodd0;
290}
291double AliD0Trigger::Energy()
292{
293 double kMD0 = 1.8645; // D0 mass
294 return sqrt(P()*P()+kMD0*kMD0);
295}
296bool AliD0Trigger::CosThetaStar()
297{
298 bool goodtheta=false;
299 double kMD0 = 1.8645; // D0 mass
300 double kMK = 0.49368; // K+ mass
301 double kMPi = 0.13957; // Pi+ mass
302 double qL=0,qL2=0,ctsD0=0,ctsD0bar=0;
303
304 Double_t pStar = TMath::Sqrt(TMath::Power(kMD0*kMD0-kMK*kMK-kMPi*kMPi,2.)-4.*kMK*kMK*kMPi*kMPi)/(2.*kMD0);
305
306 Double_t beta = P()/Energy();
307 Double_t gamma = Energy()/kMD0;
308
309 TVector3 mom(momenta[0],momenta[1],momenta[2]);
310 TVector3 mom2(momenta[3],momenta[4],momenta[5]);
311 TVector3 momD(Px(),Py(),Pz());
312
313 qL = mom.Dot(momD)/momD.Mag();
314
315 ctsD0 = (qL/gamma-beta*TMath::Sqrt(pStar*pStar+kMK*kMK))/pStar;
316
317 qL2 = mom.Dot(momD)/momD.Mag();
318
319 ctsD0bar = (qL2/gamma-beta*TMath::Sqrt(pStar*pStar+kMK*kMK))/pStar;
320
321 if(TMath::Abs(ctsD0) > cutCosThetaStar){
322 if(TMath::Abs(ctsD0bar) > cutCosThetaStar){
323 goodtheta=true;
324 }
325 }
326 return goodtheta;
327}
328bool AliD0Trigger::pTchild()
329{
330 bool goodpT=false;
331 double pT1=0,pT2=0;
332
333 pT1=sqrt(pow(momenta[0],2)+pow(momenta[1],2));
334 pT2=sqrt(pow(momenta[3],2)+pow(momenta[4],2));
335
336 if(pT1>cutpTchild){
337 if(pT2>cutpTchild){
338 goodpT=true;
339 }
340 }
341 return goodpT;
342}