[u/mrichter/AliRoot.git] / HLT / trigger / AliD0Trigger.cxx

#include "AliD0Trigger.h"
#include <TMath.h>
#include "AliITStrackV2.h"
#include "AliL3Transform.h"
#include <TVector3.h>

ClassImp(AliD0Trigger)

AliD0Trigger::AliD0Trigger()
{
  posTrack=0;
  negTrack=0;
}
AliD0Trigger::AliD0Trigger(double c[4],double b,double pv[3])
{
  posTrack=0;
  negTrack=0;
  cutV0low=c[0]; 
  cutV0high=c[1];
  cutInvMass=c[2]; 
  cutPointAngle=c[3];
  Bfield=b;
  primaryVertex[0]=pv[0];
  primaryVertex[1]=pv[1];
  primaryVertex[2]=pv[2];
}
AliD0Trigger::AliD0Trigger(AliITStrackV2 * posT, AliITStrackV2 * negT)
{ 
  posTrack=posT;
  negTrack=negT;
}


AliD0Trigger::~AliD0Trigger()
{
 
}

void AliD0Trigger::SetTracks(AliITStrackV2 * posT, AliITStrackV2 * negT){
  posTrack=posT;
  negTrack=negT;
}
bool AliD0Trigger::FindInvMass() {
   
  double mD0 = 1.8645;
  double mPi = 0.13957;
  double mKa = 0.49368;
  
  double mom2[2],momTot2,energy[2];
  
  mom2[0]=pow(momenta[0],2)+pow(momenta[1],2)+pow(momenta[2],2);
  mom2[1]=pow(momenta[3],2)+pow(momenta[4],2)+pow(momenta[5],2);
  
  momTot2 = pow(momenta[0]+momenta[3],2)+pow(momenta[1]+momenta[4],2)+pow(momenta[2]+momenta[5],2);
  
  // D0 -> K- Pi+
  energy[1] = sqrt(pow(mKa,2)+mom2[1]);
  energy[0] = sqrt(pow(mPi,2)+mom2[0]);
  
  double minvD0 = sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
   
  // D0bar -> K+ Pi-
  energy[0] = sqrt(mKa*mKa+mom2[0]);
  energy[1] = sqrt(mPi*mPi+mom2[1]);
  
  double minvD0bar = sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
  
  if(fabs(minvD0-mD0) < cutInvMass)    return true;
  if(fabs(minvD0bar-mD0) < cutInvMass) return true;
  return false;
  
}
bool AliD0Trigger::FindV0(){

  bool goodV0=false;	
  
  double Gxpos=posTrack->GetX()*cos(posTrack->GetAlpha())-posTrack->GetY()*sin(posTrack->GetAlpha());
  double Gypos=posTrack->GetX()*sin(posTrack->GetAlpha())+posTrack->GetY()*cos(posTrack->GetAlpha());
  double Gxneg=negTrack->GetX()*cos(negTrack->GetAlpha())-negTrack->GetY()*sin(negTrack->GetAlpha());
  double Gyneg=negTrack->GetX()*sin(negTrack->GetAlpha())+negTrack->GetY()*cos(negTrack->GetAlpha());
    
  double r1=fabs(1/(AliL3Transform::GetBFact()*Bfield*posTrack->Get1Pt()));
  double r2=fabs(1/(AliL3Transform::GetBFact()*Bfield*negTrack->Get1Pt()));
  //double a1=posTrack->GetX()-(r1*cos(asin(posTrack->GetSnp())+TMath::PiOver2()));
  //double a2=negTrack->GetX()-(r1*cos(asin(negTrack->GetSnp())-TMath::PiOver2()));
  //double b1=posTrack->GetY()-(r1*sin(asin(posTrack->GetSnp())+TMath::PiOver2()));
  //double b2=negTrack->GetY()-(r1*sin(asin(negTrack->GetSnp())-TMath::PiOver2()));
  
  double a1=Gxpos-(r1*cos(posTrack->GetAlpha()+asin(posTrack->GetSnp())+TMath::PiOver2()));
  double a2=Gxneg-(r2*cos(negTrack->GetAlpha()+asin(negTrack->GetSnp())-TMath::PiOver2()));
  double b1=Gypos-(r1*sin(posTrack->GetAlpha()+asin(posTrack->GetSnp())+TMath::PiOver2()));
  double b2=Gyneg-(r2*sin(negTrack->GetAlpha()+asin(negTrack->GetSnp())-TMath::PiOver2()));

  double a3=a1-a2;
  double b3=b1-b2;
  double r3=sqrt(pow(a3,2)+pow(b3,2));
  
  if(r3<r1+r2){
    double q=a1-a2;
    double p=2*(b1-b2);
    double t=pow(r2,2)-pow(r2,2)+pow(b1,2)-pow(b2,2)-pow(q,2);
 
    double Fa=pow(p,2)+(4*pow(q,2));
    double Fb=-(2*t*p+8*pow(q,2)*b1);
    double Fc=pow(t,2)-(4*pow(q,2)*pow(r1,2))+(4*pow(q,2)*pow(b1,2));
 
    if(pow(Fb,2)-(4*Fa*Fc)>=0){
      double y1=(-Fb+(sqrt(pow(Fb,2)-(4*Fa*Fc))))/(2*Fa);  //noe feil her. floating point
      double y2=(-Fb-(sqrt(pow(Fb,2)-(4*Fa*Fc))))/(2*Fa);  //trolig negativ under rot 
      
      double x1=sqrt(pow(r1,2)-pow((y1-b1),2))+a1;    
      double x2=sqrt(pow(r1,2)-pow((y2-b1),2))+a1;    
      double x3=-sqrt(pow(r1,2)-pow((y1-b1),2))+a1;
      double x4=-sqrt(pow(r1,2)-pow((y2-b1),2))+a1;
      
      double V01=sqrt(pow(x1,2)+pow(y1,2));
      double V02=sqrt(pow(x2,2)+pow(y2,2));
      double V03=sqrt(pow(x3,2)+pow(y1,2));
      double V04=sqrt(pow(x4,2)+pow(y2,2));
      
      double V0[4]={V01,V02,V03,V04};
      int index=0;
      
      double nearestV0=V0[0];
      for(int i=1;i<4;i++){
	if(nearestV0>V0[i]){
	  nearestV0=V0[i];
	  index=i;
	}
      }
      if(nearestV0<cutV0high && nearestV0>cutV0low){  
      	goodV0=true;
	switch (index){
	case 0 : bestV0[0]=x1;bestV0[1]=y1;bestV0[0]=0.;break;
	case 1 : bestV0[0]=x2;bestV0[1]=y2;bestV0[0]=0.;break;
	case 2 : bestV0[0]=x3;bestV0[1]=y1;bestV0[0]=0.;break;
	case 3 : bestV0[0]=x4;bestV0[1]=y2;bestV0[0]=0.;break;
	default : printf("Can't set V0");
	}
      }
    }
  }
  return goodV0;	
}
void AliD0Trigger::FindMomentaAtVertex(){

  double r1=fabs(1/(AliL3Transform::GetBFact()*Bfield*posTrack->Get1Pt()));
  double r2=fabs(1/(AliL3Transform::GetBFact()*Bfield*negTrack->Get1Pt()));
  
  double Gx1=posTrack->GetX()*cos(posTrack->GetAlpha())-posTrack->GetY()*sin(posTrack->GetAlpha());
  double Gy1=posTrack->GetX()*sin(posTrack->GetAlpha())+posTrack->GetY()*cos(posTrack->GetAlpha());
  double Gx2=negTrack->GetX()*cos(negTrack->GetAlpha())-negTrack->GetY()*sin(negTrack->GetAlpha());
  double Gy2=negTrack->GetX()*sin(negTrack->GetAlpha())+negTrack->GetY()*cos(negTrack->GetAlpha());
  
  double centerx1=Gx1-(r1*cos(posTrack->GetAlpha()+asin(posTrack->GetSnp())+TMath::PiOver2()));
  double centerx2=Gx2-(r2*cos(negTrack->GetAlpha()+asin(negTrack->GetSnp())-TMath::PiOver2()));
  double centery1=Gy1-(r1*sin(posTrack->GetAlpha()+asin(posTrack->GetSnp())+TMath::PiOver2()));
  double centery2=Gy2-(r2*sin(negTrack->GetAlpha()+asin(negTrack->GetSnp())-TMath::PiOver2()));
  
  double a1=sqrt(pow(Gx1-bestV0[0],2)+pow(Gy1-bestV0[1],2));
  double b1=sqrt(pow(bestV0[0]-centerx1,2)+pow(bestV0[1]-centery1,2));
  double c1=sqrt(pow(Gx1-centerx1,2)+pow(Gy1-centery1,2));
  double a2=sqrt(pow(Gx2-bestV0[0],2)+pow(Gy2-bestV0[1],2));
  double b2=sqrt(pow(bestV0[0]-centerx2,2)+pow(bestV0[1]-centery2,2));
  double c2=sqrt(pow(Gx2-centerx2,2)+pow(Gy2-centery2,2));
  
  double alpha1=acos((pow(b1,2)+pow(c1,2)-pow(a1,2))/(2*b1*c1));
  double alpha2=acos((pow(b2,2)+pow(c2,2)-pow(a2,2))/(2*b2*c2));

  double sign1=0;
  double sign2=0;
  
  if((sqrt(pow(Gx1,2)+pow(Gy1,2)))>(sqrt(pow(bestV0[0],2)+pow(bestV0[1],2))))
    sign1=-1;
  else
    sign1=1;
  if((sqrt(pow(Gx2,2)+pow(Gy2,2)))>(sqrt(pow(bestV0[0],2)+pow(bestV0[1],2))))
    sign2=1;
  else
    sign2=-1;
  
  double psi1=(posTrack->GetAlpha()+asin(posTrack->GetSnp()))+(sign1*alpha1);
  double psi2=(negTrack->GetAlpha()+asin(negTrack->GetSnp()))+(sign2*alpha2);
  
  double ptP = 1./fabs(posTrack->Get1Pt());
  momenta[0] = ptP*cos(psi1); 
  momenta[1] = ptP*sin(psi1);
  momenta[2] = ptP*posTrack->GetTgl();
  
  double ptN = 1./fabs(negTrack->Get1Pt());
  momenta[3] = ptN*cos(psi2); 
  momenta[4] = ptN*sin(psi2);
  momenta[5] = ptN*negTrack->GetTgl();
  
}
Bool_t AliD0Trigger::PointingAngle(){

  TVector3 mom(momenta[0]+momenta[3],momenta[1]+momenta[4],momenta[2]+momenta[5]);
  TVector3 flight(bestV0[0]-primaryVertex[0],bestV0[1]-primaryVertex[1],bestV0[2]-primaryVertex[2]);

  double pta = mom.Angle(flight);

  if(cos(pta)<cutPointAngle)
    return true;
  else
    return false;
}
Commit	Line	Data
5a31e9df	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()BfieldposTrack->Get1Pt()));
83	double r2=fabs(1/(AliL3Transform::GetBFact()BfieldnegTrack->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=-(2tp+8pow(q,2)b1);
105	double Fc=pow(t,2)-(4pow(q,2)pow(r1,2))+(4pow(q,2)pow(b1,2));
106
107	if(pow(Fb,2)-(4FaFc)>=0){
108	double y1=(-Fb+(sqrt(pow(Fb,2)-(4FaFc))))/(2*Fa); //noe feil her. floating point
109	double y2=(-Fb-(sqrt(pow(Fb,2)-(4FaFc))))/(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()BfieldposTrack->Get1Pt()));
148	double r2=fabs(1/(AliL3Transform::GetBFact()BfieldnegTrack->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))/(2b1c1));
168	double alpha2=acos((pow(b2,2)+pow(c2,2)-pow(a2,2))/(2b2c2));
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	}