Added support for static libraries,
[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>
6
7ClassImp(AliD0Trigger)
8
9AliD0Trigger::AliD0Trigger()
10{
11 posTrack=0;
12 negTrack=0;
13}
14AliD0Trigger::AliD0Trigger(double c[4],double b,double pv[3])
15{
16 posTrack=0;
17 negTrack=0;
18 cutV0low=c[0];
19 cutV0high=c[1];
20 cutInvMass=c[2];
21 cutPointAngle=c[3];
22 Bfield=b;
23 primaryVertex[0]=pv[0];
24 primaryVertex[1]=pv[1];
25 primaryVertex[2]=pv[2];
26}
27AliD0Trigger::AliD0Trigger(AliITStrackV2 * posT, AliITStrackV2 * negT)
28{
29 posTrack=posT;
30 negTrack=negT;
31}
32
33
34AliD0Trigger::~AliD0Trigger()
35{
36
37}
38
39void AliD0Trigger::SetTracks(AliITStrackV2 * posT, AliITStrackV2 * negT){
40 posTrack=posT;
41 negTrack=negT;
42}
43bool AliD0Trigger::FindInvMass() {
44
45 double mD0 = 1.8645;
46 double mPi = 0.13957;
47 double mKa = 0.49368;
48
49 double mom2[2],momTot2,energy[2];
50
51 mom2[0]=pow(momenta[0],2)+pow(momenta[1],2)+pow(momenta[2],2);
52 mom2[1]=pow(momenta[3],2)+pow(momenta[4],2)+pow(momenta[5],2);
53
54 momTot2 = pow(momenta[0]+momenta[3],2)+pow(momenta[1]+momenta[4],2)+pow(momenta[2]+momenta[5],2);
55
56 // D0 -> K- Pi+
57 energy[1] = sqrt(pow(mKa,2)+mom2[1]);
58 energy[0] = sqrt(pow(mPi,2)+mom2[0]);
59
60 double minvD0 = sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
61
62 // D0bar -> K+ Pi-
63 energy[0] = sqrt(mKa*mKa+mom2[0]);
64 energy[1] = sqrt(mPi*mPi+mom2[1]);
65
66 double minvD0bar = sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
67
68 if(fabs(minvD0-mD0) < cutInvMass) return true;
69 if(fabs(minvD0bar-mD0) < cutInvMass) return true;
70 return false;
71
72}
73bool AliD0Trigger::FindV0(){
74
75 bool goodV0=false;
76
77 double Gxpos=posTrack->GetX()*cos(posTrack->GetAlpha())-posTrack->GetY()*sin(posTrack->GetAlpha());
78 double Gypos=posTrack->GetX()*sin(posTrack->GetAlpha())+posTrack->GetY()*cos(posTrack->GetAlpha());
79 double Gxneg=negTrack->GetX()*cos(negTrack->GetAlpha())-negTrack->GetY()*sin(negTrack->GetAlpha());
80 double Gyneg=negTrack->GetX()*sin(negTrack->GetAlpha())+negTrack->GetY()*cos(negTrack->GetAlpha());
81
82 double r1=fabs(1/(AliL3Transform::GetBFact()*Bfield*posTrack->Get1Pt()));
83 double r2=fabs(1/(AliL3Transform::GetBFact()*Bfield*negTrack->Get1Pt()));
84 //double a1=posTrack->GetX()-(r1*cos(asin(posTrack->GetSnp())+TMath::PiOver2()));
85 //double a2=negTrack->GetX()-(r1*cos(asin(negTrack->GetSnp())-TMath::PiOver2()));
86 //double b1=posTrack->GetY()-(r1*sin(asin(posTrack->GetSnp())+TMath::PiOver2()));
87 //double b2=negTrack->GetY()-(r1*sin(asin(negTrack->GetSnp())-TMath::PiOver2()));
88
89 double a1=Gxpos-(r1*cos(posTrack->GetAlpha()+asin(posTrack->GetSnp())+TMath::PiOver2()));
90 double a2=Gxneg-(r2*cos(negTrack->GetAlpha()+asin(negTrack->GetSnp())-TMath::PiOver2()));
91 double b1=Gypos-(r1*sin(posTrack->GetAlpha()+asin(posTrack->GetSnp())+TMath::PiOver2()));
92 double b2=Gyneg-(r2*sin(negTrack->GetAlpha()+asin(negTrack->GetSnp())-TMath::PiOver2()));
93
94 double a3=a1-a2;
95 double b3=b1-b2;
96 double r3=sqrt(pow(a3,2)+pow(b3,2));
97
98 if(r3<r1+r2){
99 double q=a1-a2;
100 double p=2*(b1-b2);
101 double t=pow(r2,2)-pow(r2,2)+pow(b1,2)-pow(b2,2)-pow(q,2);
102
103 double Fa=pow(p,2)+(4*pow(q,2));
104 double Fb=-(2*t*p+8*pow(q,2)*b1);
105 double Fc=pow(t,2)-(4*pow(q,2)*pow(r1,2))+(4*pow(q,2)*pow(b1,2));
106
107 if(pow(Fb,2)-(4*Fa*Fc)>=0){
108 double y1=(-Fb+(sqrt(pow(Fb,2)-(4*Fa*Fc))))/(2*Fa); //noe feil her. floating point
109 double y2=(-Fb-(sqrt(pow(Fb,2)-(4*Fa*Fc))))/(2*Fa); //trolig negativ under rot
110
111 double x1=sqrt(pow(r1,2)-pow((y1-b1),2))+a1;
112 double x2=sqrt(pow(r1,2)-pow((y2-b1),2))+a1;
113 double x3=-sqrt(pow(r1,2)-pow((y1-b1),2))+a1;
114 double x4=-sqrt(pow(r1,2)-pow((y2-b1),2))+a1;
115
116 double V01=sqrt(pow(x1,2)+pow(y1,2));
117 double V02=sqrt(pow(x2,2)+pow(y2,2));
118 double V03=sqrt(pow(x3,2)+pow(y1,2));
119 double V04=sqrt(pow(x4,2)+pow(y2,2));
120
121 double V0[4]={V01,V02,V03,V04};
122 int index=0;
123
124 double nearestV0=V0[0];
125 for(int i=1;i<4;i++){
126 if(nearestV0>V0[i]){
127 nearestV0=V0[i];
128 index=i;
129 }
130 }
131 if(nearestV0<cutV0high && nearestV0>cutV0low){
132 goodV0=true;
133 switch (index){
134 case 0 : bestV0[0]=x1;bestV0[1]=y1;bestV0[0]=0.;break;
135 case 1 : bestV0[0]=x2;bestV0[1]=y2;bestV0[0]=0.;break;
136 case 2 : bestV0[0]=x3;bestV0[1]=y1;bestV0[0]=0.;break;
137 case 3 : bestV0[0]=x4;bestV0[1]=y2;bestV0[0]=0.;break;
138 default : printf("Can't set V0");
139 }
140 }
141 }
142 }
143 return goodV0;
144}
145void AliD0Trigger::FindMomentaAtVertex(){
146
147 double r1=fabs(1/(AliL3Transform::GetBFact()*Bfield*posTrack->Get1Pt()));
148 double r2=fabs(1/(AliL3Transform::GetBFact()*Bfield*negTrack->Get1Pt()));
149
150 double Gx1=posTrack->GetX()*cos(posTrack->GetAlpha())-posTrack->GetY()*sin(posTrack->GetAlpha());
151 double Gy1=posTrack->GetX()*sin(posTrack->GetAlpha())+posTrack->GetY()*cos(posTrack->GetAlpha());
152 double Gx2=negTrack->GetX()*cos(negTrack->GetAlpha())-negTrack->GetY()*sin(negTrack->GetAlpha());
153 double Gy2=negTrack->GetX()*sin(negTrack->GetAlpha())+negTrack->GetY()*cos(negTrack->GetAlpha());
154
155 double centerx1=Gx1-(r1*cos(posTrack->GetAlpha()+asin(posTrack->GetSnp())+TMath::PiOver2()));
156 double centerx2=Gx2-(r2*cos(negTrack->GetAlpha()+asin(negTrack->GetSnp())-TMath::PiOver2()));
157 double centery1=Gy1-(r1*sin(posTrack->GetAlpha()+asin(posTrack->GetSnp())+TMath::PiOver2()));
158 double centery2=Gy2-(r2*sin(negTrack->GetAlpha()+asin(negTrack->GetSnp())-TMath::PiOver2()));
159
160 double a1=sqrt(pow(Gx1-bestV0[0],2)+pow(Gy1-bestV0[1],2));
161 double b1=sqrt(pow(bestV0[0]-centerx1,2)+pow(bestV0[1]-centery1,2));
162 double c1=sqrt(pow(Gx1-centerx1,2)+pow(Gy1-centery1,2));
163 double a2=sqrt(pow(Gx2-bestV0[0],2)+pow(Gy2-bestV0[1],2));
164 double b2=sqrt(pow(bestV0[0]-centerx2,2)+pow(bestV0[1]-centery2,2));
165 double c2=sqrt(pow(Gx2-centerx2,2)+pow(Gy2-centery2,2));
166
167 double alpha1=acos((pow(b1,2)+pow(c1,2)-pow(a1,2))/(2*b1*c1));
168 double alpha2=acos((pow(b2,2)+pow(c2,2)-pow(a2,2))/(2*b2*c2));
169
170 double sign1=0;
171 double sign2=0;
172
173 if((sqrt(pow(Gx1,2)+pow(Gy1,2)))>(sqrt(pow(bestV0[0],2)+pow(bestV0[1],2))))
174 sign1=-1;
175 else
176 sign1=1;
177 if((sqrt(pow(Gx2,2)+pow(Gy2,2)))>(sqrt(pow(bestV0[0],2)+pow(bestV0[1],2))))
178 sign2=1;
179 else
180 sign2=-1;
181
182 double psi1=(posTrack->GetAlpha()+asin(posTrack->GetSnp()))+(sign1*alpha1);
183 double psi2=(negTrack->GetAlpha()+asin(negTrack->GetSnp()))+(sign2*alpha2);
184
185 double ptP = 1./fabs(posTrack->Get1Pt());
186 momenta[0] = ptP*cos(psi1);
187 momenta[1] = ptP*sin(psi1);
188 momenta[2] = ptP*posTrack->GetTgl();
189
190 double ptN = 1./fabs(negTrack->Get1Pt());
191 momenta[3] = ptN*cos(psi2);
192 momenta[4] = ptN*sin(psi2);
193 momenta[5] = ptN*negTrack->GetTgl();
194
195}
196Bool_t AliD0Trigger::PointingAngle(){
197
198 TVector3 mom(momenta[0]+momenta[3],momenta[1]+momenta[4],momenta[2]+momenta[5]);
199 TVector3 flight(bestV0[0]-primaryVertex[0],bestV0[1]-primaryVertex[1],bestV0[2]-primaryVertex[2]);
200
201 double pta = mom.Angle(flight);
202
203 if(cos(pta)<cutPointAngle)
204 return true;
205 else
206 return false;
207}