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

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

ClassImp(AliD0Trigger)

AliD0Trigger::AliD0Trigger()
{
  posTrack=0;
  negTrack=0;
}
AliD0Trigger::AliD0Trigger(double c[6],double b,double pv[3])
{
  posTrack=0;
  negTrack=0;
  cutV0low=c[0]; 
  cutV0high=c[1];
  cutInvMass=c[2]; 
  cutPointAngle=c[3];
  cutd0d0=c[4];
  cutCosThetaStar=c[5];
  cutpTchild=c[6];
  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=0,Gypos=0,Gxneg=0,Gyneg=0;
  double r1=0,r2=0,a1=0,a2=0,b1=0,b2=0;
  double q=0,p=0,t=0,Fa=0,Fb=0,Fc=0;
  double y1=0,y2=0,x1=0,x2=0,x3=0,x4=0;
  double V01=0,V02=0,V03=0,V04=0;
  double V0[4]={0,0,0,0};
  bestV0[0]=0;  bestV0[1]=0;  bestV0[2]=0;
  
  //Do not get right values
  Gxpos=posTrack->GetX()*cos(posTrack->GetAlpha())-posTrack->GetY()*sin(posTrack->GetAlpha());
  Gypos=posTrack->GetX()*sin(posTrack->GetAlpha())+posTrack->GetY()*cos(posTrack->GetAlpha());
  Gxneg=negTrack->GetX()*cos(negTrack->GetAlpha())-negTrack->GetY()*sin(negTrack->GetAlpha());
  Gyneg=negTrack->GetX()*sin(negTrack->GetAlpha())+negTrack->GetY()*cos(negTrack->GetAlpha());
    
  r1=fabs(1/(AliL3Transform::GetBFact()*Bfield*posTrack->Get1Pt()));
  r2=fabs(1/(AliL3Transform::GetBFact()*Bfield*negTrack->Get1Pt()));
  
  a1=Gxpos-(r1*cos(posTrack->GetAlpha()+asin(posTrack->GetSnp())+TMath::PiOver2()));
  a2=Gxneg-(r2*cos(negTrack->GetAlpha()+asin(negTrack->GetSnp())-TMath::PiOver2()));
  b1=Gypos-(r1*sin(posTrack->GetAlpha()+asin(posTrack->GetSnp())+TMath::PiOver2()));
  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){
  q=a1-a2;
  p=2*(b1-b2);
  t=pow(r2,2)-pow(r2,2)+pow(b1,2)-pow(b2,2)-pow(q,2);
  
  Fa=pow(p,2)+(4*pow(q,2));
  Fb=-(2*t*p+8*pow(q,2)*b1);
  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){
    y1=(-Fb+(sqrt(pow(Fb,2)-(4*Fa*Fc))))/(2*Fa);  
    y2=(-Fb-(sqrt(pow(Fb,2)-(4*Fa*Fc))))/(2*Fa);  
    
    x1=sqrt(pow(r1,2)-pow((y1-b1),2))+a1;    
    x2=sqrt(pow(r1,2)-pow((y2-b1),2))+a1;    
    x3=-sqrt(pow(r1,2)-pow((y1-b1),2))+a1;
    x4=-sqrt(pow(r1,2)-pow((y2-b1),2))+a1;
    
    V01=sqrt(pow(x1-primaryVertex[0],2)+pow(y1-primaryVertex[1],2));
    V02=sqrt(pow(x2-primaryVertex[0],2)+pow(y2-primaryVertex[1],2));
    V03=sqrt(pow(x3-primaryVertex[0],2)+pow(y1-primaryVertex[1],2));  
    V04=sqrt(pow(x4-primaryVertex[0],2)+pow(y2-primaryVertex[1],2));  
    
    V0[0]=V01;V0[1]=V02;V0[2]=V03;V0[3]=V04;

    int index=0;
    
    double nearestV0=V0[0];

    for(int i=0;i<4;i++){                        
      if(nearestV0>V0[i]){
	nearestV0=V0[i];
	index=i;
      }
    }
    cout<<"index: "<<index<<endl;
    if(nearestV0<cutV0high){
      if(nearestV0>cutV0low){  
	goodV0=true;
	if(index==0){
	  bestV0[0]=x1;   bestV0[1]=y1;   bestV0[2]=0.;
	}
	if(index==1){
	  bestV0[0]=x2;   bestV0[1]=y2;   bestV0[2]=0.;
	}
	if(index==2){
	  bestV0[0]=x3;   bestV0[1]=y1;   bestV0[2]=0.;
	}
	if(index==3){
	  bestV0[0]=x4;   bestV0[1]=y2;   bestV0[2]=0.;
	}
      }
    }
    //}
  }
  
  cout<<"My V0: x: "<<bestV0[0]<<" y: "<<bestV0[1]<<" z: "<<bestV0[2]<<endl;

  return goodV0;	
}
void AliD0Trigger::FindMomentaAtVertex(){

  //This method moves the momenta to the secondary vertex

  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(Px(),Py(),Pz());
  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;
}
void AliD0Trigger::SetMomenta(double m[6])
{
  momenta[0] = m[0]; 
  momenta[1] = m[1]; 
  momenta[2] = m[2]; 
  momenta[3] = m[3]; 
  momenta[4] = m[4]; 
  momenta[5] = m[5]; 
}
void AliD0Trigger::FindMomentaOffline()
{
  double ptP,alphaP,phiP,ptN,alphaN,phiN;
  // momenta of the tracks at the vertex
  ptP = 1./TMath::Abs(posTrack->Get1Pt());
  alphaP = posTrack->GetAlpha();
  phiP = alphaP+TMath::ASin(posTrack->GetSnp());
  momenta[0] = ptP*TMath::Cos(phiP); 
  momenta[1] = ptP*TMath::Sin(phiP);
  momenta[2] = ptP*posTrack->GetTgl();
  
  ptN = 1./TMath::Abs(negTrack->Get1Pt());
  alphaN = negTrack->GetAlpha();
  phiN = alphaN+TMath::ASin(negTrack->GetSnp());
  momenta[3] = ptN*TMath::Cos(phiN); 
  momenta[4] = ptN*TMath::Sin(phiN);
  momenta[5] = ptN*negTrack->GetTgl();  
}
bool AliD0Trigger::FindV0offline(double v[3])
{

  bool goodV0=false;
  double r=0;
  
  bestV0[0]=v[0];
  bestV0[1]=v[1];
  bestV0[2]=v[2];

  r=sqrt(pow(bestV0[0]-primaryVertex[0],2)+pow(bestV0[1]-primaryVertex[1],2)+pow(bestV0[2]-primaryVertex[2],2));

  if(r<cutV0high){
    if(r>cutV0low){
      goodV0=true;
    }
  }

  return goodV0;

}
bool AliD0Trigger::d0d0()
{
  bool goodd0=false;
  double d00=0, d01=0;

  d00 =  10000.*posTrack->GetD(primaryVertex[0],primaryVertex[1]);
  d01 = -10000.*negTrack->GetD(primaryVertex[0],primaryVertex[1]);

  if(d00*d01<cutd0d0){
    goodd0=true;
  }

  return goodd0;
}
double AliD0Trigger::Energy()
{
  double kMD0 = 1.8645;  // D0  mass
  return sqrt(P()*P()+kMD0*kMD0);
}
bool AliD0Trigger::CosThetaStar()
{
  bool goodtheta=false;
  double kMD0 = 1.8645;  // D0  mass
  double kMK  = 0.49368; // K+  mass
  double kMPi = 0.13957; // Pi+ mass
  double qL=0,qL2=0,ctsD0=0,ctsD0bar=0;

  Double_t pStar = TMath::Sqrt(TMath::Power(kMD0*kMD0-kMK*kMK-kMPi*kMPi,2.)-4.*kMK*kMK*kMPi*kMPi)/(2.*kMD0);

  Double_t beta = P()/Energy();
  Double_t gamma = Energy()/kMD0;

  TVector3 mom(momenta[0],momenta[1],momenta[2]);
  TVector3 mom2(momenta[3],momenta[4],momenta[5]);
  TVector3 momD(Px(),Py(),Pz());

  qL = mom.Dot(momD)/momD.Mag();

  ctsD0 = (qL/gamma-beta*TMath::Sqrt(pStar*pStar+kMK*kMK))/pStar;

  qL2 = mom.Dot(momD)/momD.Mag();

  ctsD0bar = (qL2/gamma-beta*TMath::Sqrt(pStar*pStar+kMK*kMK))/pStar;

  if(fabs(ctsD0) > cutCosThetaStar || fabs(ctsD0bar) > cutCosThetaStar){
      goodtheta=true;
    }
  return goodtheta;
}
bool AliD0Trigger::pTchild()
{
  bool goodpT=false;
  double pT1=0,pT2=0;
  
  pT1=sqrt(pow(momenta[0],2)+pow(momenta[1],2));
  pT2=sqrt(pow(momenta[3],2)+pow(momenta[4],2));

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