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