Added support for static libraries,
[u/mrichter/AliRoot.git] / HLT / trigger / AliD0Trigger.cxx
1 #include "AliD0Trigger.h"
2 #include <TMath.h>
3 #include "AliITStrackV2.h"
4 #include "AliL3Transform.h"
5 #include <TVector3.h>
6
7 ClassImp(AliD0Trigger)
8
9 AliD0Trigger::AliD0Trigger()
10 {
11   posTrack=0;
12   negTrack=0;
13 }
14 AliD0Trigger::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 }
27 AliD0Trigger::AliD0Trigger(AliITStrackV2 * posT, AliITStrackV2 * negT)
28
29   posTrack=posT;
30   negTrack=negT;
31 }
32
33
34 AliD0Trigger::~AliD0Trigger()
35 {
36  
37 }
38
39 void AliD0Trigger::SetTracks(AliITStrackV2 * posT, AliITStrackV2 * negT){
40   posTrack=posT;
41   negTrack=negT;
42 }
43 bool 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 }
73 bool 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 }
145 void 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 }
196 Bool_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 }