[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
	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	}