Added Gautes changes from Bergen.
[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 #include <iostream.h>
7
8 ClassImp(AliD0Trigger)
9
10 AliD0Trigger::AliD0Trigger()
11 {
12   posTrack=0;
13   negTrack=0;
14 }
15 AliD0Trigger::AliD0Trigger(double c[7],double b,double pv[3])
16 {
17   posTrack=0;
18   negTrack=0;
19   cutV0low=c[0]; 
20   cutV0high=c[1];
21   cutInvMass=c[2]; 
22   cutPointAngle=c[3];
23   cutd0d0=c[4];
24   cutCosThetaStar=c[5];
25   cutpTchild=c[6];
26   Bfield=b;
27   primaryVertex[0]=pv[0];
28   primaryVertex[1]=pv[1];
29   primaryVertex[2]=pv[2];
30 }
31 AliD0Trigger::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 }
47 AliD0Trigger::AliD0Trigger(AliITStrackV2 * posT, AliITStrackV2 * negT)
48
49   posTrack=posT;
50   negTrack=negT;
51 }
52
53
54 AliD0Trigger::~AliD0Trigger()
55 {
56  
57 }
58
59 void AliD0Trigger::SetTracks(AliITStrackV2 * posT, AliITStrackV2 * negT){
60   posTrack=posT;
61   negTrack=negT;
62 }
63 bool 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 }
93 bool AliD0Trigger::FindV0(){
94
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;
104   
105   //Do not get right values
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());
110     
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){
133     y1=(-Fb+(sqrt(pow(Fb,2)-(4*Fa*Fc))))/(2*Fa);  
134     y2=(-Fb-(sqrt(pow(Fb,2)-(4*Fa*Fc))))/(2*Fa);  
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;
156       }
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.;
173         }
174       }
175     }
176     //}
177   }
178   
179   cout<<"My V0: x: "<<bestV0[0]<<" y: "<<bestV0[1]<<" z: "<<bestV0[2]<<endl;
180
181   return goodV0;        
182 }
183 void AliD0Trigger::FindMomentaAtVertex(){
184
185   //This method moves the momenta to the secondary vertex
186
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 }
236 Bool_t AliD0Trigger::PointingAngle(){
237
238   TVector3 mom(Px(),Py(),Pz());
239   TVector3 flight(bestV0[0]-primaryVertex[0],bestV0[1]-primaryVertex[1],bestV0[2]-primaryVertex[2]);
240
241   double pta = mom.Angle(flight);
242
243   if(cos(pta)>cutPointAngle)
244     return true;
245   else
246     return false;
247 }
248 void 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 }
257 void 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 }
275 bool 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 }
296 bool 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 }
310 double AliD0Trigger::Energy()
311 {
312   double kMD0 = 1.8645;  // D0  mass
313   return sqrt(P()*P()+kMD0*kMD0);
314 }
315 bool 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
340   if(fabs(ctsD0) > cutCosThetaStar || fabs(ctsD0bar) > cutCosThetaStar){
341       goodtheta=true;
342     }
343   return goodtheta;
344 }
345 bool 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 }
360 void AliD0Trigger::SetV0(double v[3])
361 {
362   bestV0[0]=v[0];
363   bestV0[1]=v[1];
364   bestV0[2]=v[2];
365 }