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[6],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(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
  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;
  88
  89   Gxpos=posTrack->GetX()*cos(posTrack->GetAlpha())-posTrack->GetY()*sin(posTrack->GetAlpha());
  90   Gypos=posTrack->GetX()*sin(posTrack->GetAlpha())+posTrack->GetY()*cos(posTrack->GetAlpha());
  91   Gxneg=negTrack->GetX()*cos(negTrack->GetAlpha())-negTrack->GetY()*sin(negTrack->GetAlpha());
  92   Gyneg=negTrack->GetX()*sin(negTrack->GetAlpha())+negTrack->GetY()*cos(negTrack->GetAlpha());
  93
  94   r1=fabs(1/(AliL3Transform::GetBFact()*Bfield*posTrack->Get1Pt()));
  95   r2=fabs(1/(AliL3Transform::GetBFact()*Bfield*negTrack->Get1Pt()));
  96
  97   a1=Gxpos-(r1*cos(posTrack->GetAlpha()+asin(posTrack->GetSnp())+TMath::PiOver2()));
  98   a2=Gxneg-(r2*cos(negTrack->GetAlpha()+asin(negTrack->GetSnp())-TMath::PiOver2()));
  99   b1=Gypos-(r1*sin(posTrack->GetAlpha()+asin(posTrack->GetSnp())+TMath::PiOver2()));
 100   b2=Gyneg-(r2*sin(negTrack->GetAlpha()+asin(negTrack->GetSnp())-TMath::PiOver2()));
 101
 102   //double a3=a1-a2;
 103   //double b3=b1-b2;
 104   //double r3=sqrt(pow(a3,2)+pow(b3,2));
 105
 106   //if(r3<r1+r2){
 107   q=a1-a2;
 108   p=2*(b1-b2);
 109   t=pow(r2,2)-pow(r2,2)+pow(b1,2)-pow(b2,2)-pow(q,2);
 110
 111   Fa=pow(p,2)+(4*pow(q,2));
 112   Fb=-(2*t*p+8*pow(q,2)*b1);
 113   Fc=pow(t,2)-(4*pow(q,2)*pow(r1,2))+(4*pow(q,2)*pow(b1,2));
 114
 115   if(pow(Fb,2)-(4*Fa*Fc)>=0){
 116     y1=(-Fb+(sqrt(pow(Fb,2)-(4*Fa*Fc))))/(2*Fa);  //noe feil her. floating point
 117     y2=(-Fb-(sqrt(pow(Fb,2)-(4*Fa*Fc))))/(2*Fa);  //trolig negativ under rot
 118
 119     x1=sqrt(pow(r1,2)-pow((y1-b1),2))+a1;
 120     x2=sqrt(pow(r1,2)-pow((y2-b1),2))+a1;
 121     x3=-sqrt(pow(r1,2)-pow((y1-b1),2))+a1;
 122     x4=-sqrt(pow(r1,2)-pow((y2-b1),2))+a1;
 123
 124     V01=sqrt(pow(x1-primaryVertex[0],2)+pow(y1-primaryVertex[1],2));
 125     V02=sqrt(pow(x2-primaryVertex[0],2)+pow(y2-primaryVertex[1],2));
 126     V03=sqrt(pow(x3-primaryVertex[0],2)+pow(y1-primaryVertex[1],2));
 127     V04=sqrt(pow(x4-primaryVertex[0],2)+pow(y2-primaryVertex[1],2));
 128
 129     V0[0]=V01;V0[1]=V02;V0[2]=V03;V0[3]=V04;
 130
 131     int index=0;
 132
 133     double nearestV0=V0[0];
 134
 135     for(int i=0;i<4;i++){
 136       if(nearestV0>V0[i]){
 137         nearestV0=V0[i];
 138         index=i;
 139       }
 140     }
 141     cout<<"index: "<<index<<endl;
 142     if(nearestV0<cutV0high){
 143       if(nearestV0>cutV0low){
 144         goodV0=true;
 145         if(index==0){
 146           bestV0[0]=x1;   bestV0[1]=y1;   bestV0[2]=0.;
 147         }
 148         if(index==1){
 149           bestV0[0]=x2;   bestV0[1]=y2;   bestV0[2]=0.;
 150         }
 151         if(index==2){
 152           bestV0[0]=x3;   bestV0[1]=y1;   bestV0[2]=0.;
 153         }
 154         if(index==3){
 155           bestV0[0]=x4;   bestV0[1]=y2;   bestV0[2]=0.;
 156         }
 157       }
 158     }
 159     //}
 160   }
 161
 162   cout<<"My V0: x: "<<bestV0[0]<<" y: "<<bestV0[1]<<" z: "<<bestV0[2]<<endl;
 163
 164   return goodV0;
 165 }
 166 void AliD0Trigger::FindMomentaAtVertex(){
 167
 168   double r1=fabs(1/(AliL3Transform::GetBFact()*Bfield*posTrack->Get1Pt()));
 169   double r2=fabs(1/(AliL3Transform::GetBFact()*Bfield*negTrack->Get1Pt()));
 170
 171   double Gx1=posTrack->GetX()*cos(posTrack->GetAlpha())-posTrack->GetY()*sin(posTrack->GetAlpha());
 172   double Gy1=posTrack->GetX()*sin(posTrack->GetAlpha())+posTrack->GetY()*cos(posTrack->GetAlpha());
 173   double Gx2=negTrack->GetX()*cos(negTrack->GetAlpha())-negTrack->GetY()*sin(negTrack->GetAlpha());
 174   double Gy2=negTrack->GetX()*sin(negTrack->GetAlpha())+negTrack->GetY()*cos(negTrack->GetAlpha());
 175
 176   double centerx1=Gx1-(r1*cos(posTrack->GetAlpha()+asin(posTrack->GetSnp())+TMath::PiOver2()));
 177   double centerx2=Gx2-(r2*cos(negTrack->GetAlpha()+asin(negTrack->GetSnp())-TMath::PiOver2()));
 178   double centery1=Gy1-(r1*sin(posTrack->GetAlpha()+asin(posTrack->GetSnp())+TMath::PiOver2()));
 179   double centery2=Gy2-(r2*sin(negTrack->GetAlpha()+asin(negTrack->GetSnp())-TMath::PiOver2()));
 180
 181   double a1=sqrt(pow(Gx1-bestV0[0],2)+pow(Gy1-bestV0[1],2));
 182   double b1=sqrt(pow(bestV0[0]-centerx1,2)+pow(bestV0[1]-centery1,2));
 183   double c1=sqrt(pow(Gx1-centerx1,2)+pow(Gy1-centery1,2));
 184   double a2=sqrt(pow(Gx2-bestV0[0],2)+pow(Gy2-bestV0[1],2));
 185   double b2=sqrt(pow(bestV0[0]-centerx2,2)+pow(bestV0[1]-centery2,2));
 186   double c2=sqrt(pow(Gx2-centerx2,2)+pow(Gy2-centery2,2));
 187
 188   double alpha1=acos((pow(b1,2)+pow(c1,2)-pow(a1,2))/(2*b1*c1));
 189   double alpha2=acos((pow(b2,2)+pow(c2,2)-pow(a2,2))/(2*b2*c2));
 190
 191   double sign1=0;
 192   double sign2=0;
 193
 194   if((sqrt(pow(Gx1,2)+pow(Gy1,2)))>(sqrt(pow(bestV0[0],2)+pow(bestV0[1],2))))
 195     sign1=-1;
 196   else
 197     sign1=1;
 198   if((sqrt(pow(Gx2,2)+pow(Gy2,2)))>(sqrt(pow(bestV0[0],2)+pow(bestV0[1],2))))
 199     sign2=1;
 200   else
 201     sign2=-1;
 202
 203   double psi1=(posTrack->GetAlpha()+asin(posTrack->GetSnp()))+(sign1*alpha1);
 204   double psi2=(negTrack->GetAlpha()+asin(negTrack->GetSnp()))+(sign2*alpha2);
 205
 206   double ptP = 1./fabs(posTrack->Get1Pt());
 207   momenta[0] = ptP*cos(psi1);
 208   momenta[1] = ptP*sin(psi1);
 209   momenta[2] = ptP*posTrack->GetTgl();
 210
 211   double ptN = 1./fabs(negTrack->Get1Pt());
 212   momenta[3] = ptN*cos(psi2);
 213   momenta[4] = ptN*sin(psi2);
 214   momenta[5] = ptN*negTrack->GetTgl();
 215
 216 }
 217 Bool_t AliD0Trigger::PointingAngle(){
 218
 219   TVector3 mom(momenta[0]+momenta[3],momenta[1]+momenta[4],momenta[2]+momenta[5]);
 220   TVector3 flight(bestV0[0]-primaryVertex[0],bestV0[1]-primaryVertex[1],bestV0[2]-primaryVertex[2]);
 221
 222   double pta = mom.Angle(flight);
 223
 224   if(cos(pta)<cutPointAngle)
 225     return true;
 226   else
 227     return false;
 228 }
 229 void AliD0Trigger::SetMomenta(double m[6])
 230 {
 231   momenta[0] = m[0];
 232   momenta[1] = m[1];
 233   momenta[2] = m[2];
 234   momenta[3] = m[3];
 235   momenta[4] = m[4];
 236   momenta[5] = m[5];
 237 }
 238 void AliD0Trigger::FindMomentaOffline()
 239 {
 240   double ptP,alphaP,phiP,ptN,alphaN,phiN;
 241   // momenta of the tracks at the vertex
 242   ptP = 1./TMath::Abs(posTrack->Get1Pt());
 243   alphaP = posTrack->GetAlpha();
 244   phiP = alphaP+TMath::ASin(posTrack->GetSnp());
 245   momenta[0] = ptP*TMath::Cos(phiP);
 246   momenta[1] = ptP*TMath::Sin(phiP);
 247   momenta[2] = ptP*posTrack->GetTgl();
 248
 249   ptN = 1./TMath::Abs(negTrack->Get1Pt());
 250   alphaN = negTrack->GetAlpha();
 251   phiN = alphaN+TMath::ASin(negTrack->GetSnp());
 252   momenta[3] = ptN*TMath::Cos(phiN);
 253   momenta[4] = ptN*TMath::Sin(phiN);
 254   momenta[5] = ptN*negTrack->GetTgl();
 255 }
 256 bool AliD0Trigger::FindV0offline(double v[3])
 257 {
 258
 259   bool goodV0=false;
 260   double r=0;
 261
 262   bestV0[0]=v[0];
 263   bestV0[1]=v[1];
 264   bestV0[2]=v[2];
 265
 266   r=sqrt(pow(bestV0[0]-primaryVertex[0],2)+pow(bestV0[1]-primaryVertex[1],2)+pow(bestV0[2]-primaryVertex[2],2));
 267
 268   if(r<cutV0high){
 269     if(r>cutV0low){
 270       goodV0=true;
 271     }
 272   }
 273
 274   return goodV0;
 275
 276 }
 277 bool AliD0Trigger::d0d0()
 278 {
 279   bool goodd0=false;
 280   double d00=0, d01=0;
 281
 282   d00 =  10000.*posTrack->GetD(primaryVertex[0],primaryVertex[1]);
 283   d01 = -10000.*negTrack->GetD(primaryVertex[0],primaryVertex[1]);
 284
 285   if(d00*d01<cutd0d0){
 286     goodd0=true;
 287   }
 288
 289   return goodd0;
 290 }
 291 double AliD0Trigger::Energy()
 292 {
 293   double kMD0 = 1.8645;  // D0  mass
 294   return sqrt(P()*P()+kMD0*kMD0);
 295 }
 296 bool AliD0Trigger::CosThetaStar()
 297 {
 298   bool goodtheta=false;
 299   double kMD0 = 1.8645;  // D0  mass
 300   double kMK  = 0.49368; // K+  mass
 301   double kMPi = 0.13957; // Pi+ mass
 302   double qL=0,qL2=0,ctsD0=0,ctsD0bar=0;
 303
 304   Double_t pStar = TMath::Sqrt(TMath::Power(kMD0*kMD0-kMK*kMK-kMPi*kMPi,2.)-4.*kMK*kMK*kMPi*kMPi)/(2.*kMD0);
 305
 306   Double_t beta = P()/Energy();
 307   Double_t gamma = Energy()/kMD0;
 308
 309   TVector3 mom(momenta[0],momenta[1],momenta[2]);
 310   TVector3 mom2(momenta[3],momenta[4],momenta[5]);
 311   TVector3 momD(Px(),Py(),Pz());
 312
 313   qL = mom.Dot(momD)/momD.Mag();
 314
 315   ctsD0 = (qL/gamma-beta*TMath::Sqrt(pStar*pStar+kMK*kMK))/pStar;
 316
 317   qL2 = mom.Dot(momD)/momD.Mag();
 318
 319   ctsD0bar = (qL2/gamma-beta*TMath::Sqrt(pStar*pStar+kMK*kMK))/pStar;
 320
 321   if(TMath::Abs(ctsD0) > cutCosThetaStar){
 322     if(TMath::Abs(ctsD0bar) > cutCosThetaStar){
 323       goodtheta=true;
 324     }
 325   }
 326   return goodtheta;
 327 }
 328 bool AliD0Trigger::pTchild()
 329 {
 330   bool goodpT=false;
 331   double pT1=0,pT2=0;
 332
 333   pT1=sqrt(pow(momenta[0],2)+pow(momenta[1],2));
 334   pT2=sqrt(pow(momenta[3],2)+pow(momenta[4],2));
 335
 336   if(pT1>cutpTchild){
 337     if(pT2>cutpTchild){
 338       goodpT=true;
 339     }
 340   }
 341   return goodpT;
 342 }