PWGHF/vertexingHF/AliAODRecoDecayHF4Prong.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 /////////////////////////////////////////////////////////////
  19 //
  20 // Base class for AOD reconstructed heavy-flavour 4-prong decay
  21 //
  22 // Authors: G.E.Bruno Giuseppe.Bruno@to.infn.it, R.Romita Rossella.Romita@ba.infn.it
  23 /////////////////////////////////////////////////////////////
  24
  25 #include <TDatabasePDG.h>
  26 #include "AliAODRecoDecayHF.h"
  27 #include "AliAODRecoDecayHF4Prong.h"
  28
  29 ClassImp(AliAODRecoDecayHF4Prong)
  30
  31 //--------------------------------------------------------------------------
  32 AliAODRecoDecayHF4Prong::AliAODRecoDecayHF4Prong() :
  33   AliAODRecoDecayHF(),
  34   fDist12toPrim(0),
  35   fDist3toPrim(0),
  36   fDist4toPrim(0)
  37 {
  38   //
  39   // Default Constructor
  40   //
  41 }
  42 //--------------------------------------------------------------------------
  43 AliAODRecoDecayHF4Prong::AliAODRecoDecayHF4Prong(AliAODVertex *vtx2,
  44                                                  Double_t *px,Double_t *py,Double_t *pz,
  45                                                  Double_t *d0,Double_t *d0err,
  46                                                  Double_t *dca, //Double_t sigvert,
  47                                                  Double_t dist12,Double_t dist3,
  48                                                  Double_t dist4,
  49                                                  Short_t charge) :
  50   AliAODRecoDecayHF(vtx2,4,charge,px,py,pz,d0,d0err),
  51   // fSigmaVert(sigvert),
  52   fDist12toPrim(dist12),
  53   fDist3toPrim(dist3),
  54   fDist4toPrim(dist4)
  55 {
  56   //
  57   // Constructor with AliAODVertex for decay vertex
  58   //
  59   SetDCAs(6,dca);
  60 }
  61 //--------------------------------------------------------------------------
  62 AliAODRecoDecayHF4Prong::AliAODRecoDecayHF4Prong(AliAODVertex *vtx2,
  63                                                  Double_t *d0,Double_t *d0err,
  64                                                  Double_t *dca, //Double_t sigvert,
  65                                                  Double_t dist12,Double_t dist3,
  66                                                  Double_t dist4,
  67                                                  Short_t charge) :
  68   AliAODRecoDecayHF(vtx2,4,charge,d0,d0err),
  69   //fSigmaVert(sigvert),
  70   fDist12toPrim(dist12),
  71   fDist3toPrim(dist3),
  72   fDist4toPrim(dist4)
  73 {
  74   //
  75   // Constructor with AliAODVertex for decay vertex and without prongs momenta
  76   //
  77   SetDCAs(6,dca);
  78 }
  79 //--------------------------------------------------------------------------
  80 AliAODRecoDecayHF4Prong::AliAODRecoDecayHF4Prong(const AliAODRecoDecayHF4Prong &source) :
  81   AliAODRecoDecayHF(source),
  82   //fSigmaVert(source.fSigmaVert),
  83   fDist12toPrim(source.fDist12toPrim),
  84   fDist3toPrim(source.fDist3toPrim),
  85   fDist4toPrim(source.fDist4toPrim)
  86 {
  87   //
  88   // Copy constructor
  89   //
  90 }
  91 //--------------------------------------------------------------------------
  92 AliAODRecoDecayHF4Prong &AliAODRecoDecayHF4Prong::operator=(const AliAODRecoDecayHF4Prong &source)
  93 {
  94   //
  95   // assignment operator
  96   //
  97   if(&source == this) return *this;
  98
  99   AliAODRecoDecayHF::operator=(source);
 100
 101   fDist12toPrim= source.fDist12toPrim;
 102   fDist3toPrim= source.fDist3toPrim;
 103   fDist4toPrim= source.fDist4toPrim;
 104   //fSigmaVert= source.fSigmaVert;
 105
 106   return *this;
 107 }
 108 //--------------------------------------------------------------------------
 109 void AliAODRecoDecayHF4Prong::InvMassD0(Double_t mD0[2]) const {
 110   //
 111   // Mass for the two D0 hypotheses
 112   //
 113   UInt_t pdg[4];
 114   pdg[0]=211; pdg[1]=321; pdg[2]=211; pdg[3]=211;
 115   mD0[0]=InvMass(4,pdg);
 116   pdg[1]=211; pdg[3]=321;
 117   mD0[1]=InvMass(4,pdg);
 118
 119   return;
 120 }
 121 //--------------------------------------------------------------------------
 122 void AliAODRecoDecayHF4Prong::InvMassD0bar(Double_t mD0bar[2]) const {
 123   //
 124   // Mass for the two D0bar hypotheses
 125   //
 126   UInt_t pdg[4];
 127   pdg[0]=321; pdg[1]=211; pdg[2]=211; pdg[3]=211;
 128   mD0bar[0]=InvMass(4,pdg);
 129   pdg[0]=211; pdg[2]=321;
 130   mD0bar[1]=InvMass(4,pdg);
 131
 132   return;
 133 }
 134 //--------------------------------------------------------------------------
 135
 136 Bool_t AliAODRecoDecayHF4Prong::SelectD0(const Double_t *cuts,Int_t &okD0,Int_t &okD0bar) const
 137 {
 138   //
 139   // This function compares the D0 with a set of cuts:
 140   //
 141   // cuts[0] = D0 invariant mass
 142   // cuts[1] = DCA between opposite sign tracks
 143   // cuts[2] = Distance between primary and two tracks vertex fDist12toPrim
 144   // cuts[3] = Distance between primary and three tracks vertex fDist3toPrim
 145   // cuts[4] = Distance between primary and two tracks vertex fDist4toPrim
 146   // cuts[5] = Cosinus of the pointing angle
 147   // cuts[6] = Transverse momentum of the D0 candidate
 148   // cuts[7] = Mass Pi+Pi- = mass of the rho0
 149   // cuts[8] = PID cut (one K in the quadruplet)
 150   //
 151   // If candidate D0 does not pass the cuts return kFALSE
 152   //
 153
 154   okD0=0; okD0bar=0;
 155   Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 156   Double_t mD0[2];
 157   Double_t mD0bar[2];
 158   InvMassD0(mD0);
 159   InvMassD0bar(mD0bar);
 160   Bool_t goodMass=kFALSE;
 161   if(TMath::Abs(mD0[0]-mD0PDG)<=cuts[0]) {goodMass=kTRUE; okD0=1;}
 162   if(TMath::Abs(mD0[1]-mD0PDG)<=cuts[0]) {goodMass=kTRUE; okD0=1;}
 163   if(TMath::Abs(mD0bar[0]-mD0PDG)<=cuts[0]) {goodMass=kTRUE; okD0bar=1;}
 164   if(TMath::Abs(mD0bar[1]-mD0PDG)<=cuts[0]) {goodMass=kTRUE; okD0bar=1;}
 165   if(!goodMass) return kFALSE;
 166
 167   //DCA btw opposite sign tracks
 168   if(cuts[1]>0.){
 169     if(GetDCA(0)>cuts[1]) return kFALSE;
 170     if(GetDCA(1)>cuts[1]) return kFALSE;
 171     if(GetDCA(2)>cuts[1]) return kFALSE;
 172     if(GetDCA(3)>cuts[1]) return kFALSE;
 173     if(GetDCA(4)>cuts[1]) return kFALSE;
 174     if(GetDCA(5)>cuts[1]) return kFALSE;
 175   }
 176   //2track cuts
 177   if(cuts[2]>0.){
 178     if(fDist12toPrim>10.)return kFALSE;
 179     if(fDist12toPrim<cuts[2])return kFALSE;
 180   }
 181   //3track cuts
 182   if(cuts[3]>0.){
 183     if(fDist3toPrim<cuts[3])return kFALSE;
 184   }
 185   //4track cuts
 186   if(cuts[4]>0.){
 187     if(fDist4toPrim<cuts[4])return kFALSE;
 188   }
 189   if(cuts[5]>-1.1){
 190     if(CosPointingAngle()<cuts[5])return kFALSE;
 191   }
 192   if(cuts[6]>0.){
 193     if(Pt()<cuts[6])return kFALSE;
 194   }
 195   if(cuts[7]>0.){
 196     Double_t massD0[2];
 197     Double_t massD0bar[2];
 198     Bool_t good=CutRhoMass(massD0,massD0bar,cuts[0],cuts[7]);
 199     if(!good) return kFALSE;
 200   }
 201
 202   return kTRUE;
 203 }
 204 //----------------------------------------------------------------------------
 205 Bool_t AliAODRecoDecayHF4Prong::CutRhoMass(Double_t massD0[2],Double_t massD0bar[2],Double_t cutMass,Double_t cutRho) const
 206 {
 207   //
 208   // Cut on rho->pipi mass for any of the pairs
 209   //
 210   Bool_t isGood=kFALSE;
 211   Int_t nprongs=4;
 212   for(Int_t i=0;i<2;i++){massD0[i]=0.;massD0bar[i]=0.;}
 213   Bool_t isRho=kFALSE;
 214   Bool_t isTrue=kFALSE;
 215   Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 216   Double_t mPDGRho=TDatabasePDG::Instance()->GetParticle(113)->Mass();
 217   Double_t minv01=InvMassRho(0,1);
 218   if(TMath::Abs(minv01-mPDGRho)<cutRho) isRho=kTRUE;
 219   if(isRho){
 220     UInt_t pdg1[4]={211,211,321,211};
 221     Double_t mass1=InvMass(nprongs,pdg1);
 222     if(TMath::Abs(mass1-mPDG)<cutMass) {isTrue=kTRUE;isGood=kTRUE;}
 223     if(isTrue) massD0bar[1]=mass1;
 224     isTrue=kFALSE;
 225     UInt_t pdg2[4]={211,211,211,321};
 226     Double_t mass2=InvMass(4,pdg2);
 227     if(TMath::Abs(mass2-mPDG)<cutMass) {isTrue=kTRUE;isGood=kTRUE;}
 228     if(isTrue) massD0[1]=mass2;
 229     isTrue=kFALSE;
 230   }
 231   Double_t minv03=InvMassRho(0,3);
 232   if(TMath::Abs(minv03-mPDGRho)<cutRho) isRho=kTRUE;
 233   if(isRho){
 234     UInt_t pdg1[4]={211,211,321,211};
 235     Double_t mass1=InvMass(4,pdg1);
 236     if(TMath::Abs(mass1-mPDG)<cutMass) {isTrue=kTRUE;isGood=kTRUE;}
 237     if(isTrue) massD0bar[1]=mass1;
 238     isTrue=kFALSE;
 239     UInt_t pdg2[4]={211,321,211,211};
 240     Double_t mass2=InvMass(4,pdg2);
 241     if(TMath::Abs(mass2-mPDG)<cutMass) {isTrue=kTRUE;isGood=kTRUE;}
 242     if(isTrue) massD0[0]=mass2;
 243     isTrue=kFALSE;
 244   }
 245   Double_t minv12=InvMassRho(1,2);
 246   if(TMath::Abs(minv12-mPDGRho)<cutRho) isRho=kTRUE;
 247   if(isRho){
 248     UInt_t pdg1[4]={321,211,211,211};
 249     Double_t mass1=InvMass(4,pdg1);
 250     if(TMath::Abs(mass1-mPDG)<cutMass) {isTrue=kTRUE;isGood=kTRUE;}
 251     if(isTrue) massD0bar[0]=mass1;
 252     isTrue=kFALSE;
 253     UInt_t pdg2[4]={211,211,211,321};
 254     Double_t mass2=InvMass(4,pdg2);
 255     if(TMath::Abs(mass2-mPDG)<cutMass) {isTrue=kTRUE;isGood=kTRUE;}
 256     if(isTrue) massD0[1]=mass2;
 257     isTrue=kFALSE;
 258   }
 259   Double_t minv23=InvMassRho(2,3);
 260   if(TMath::Abs(minv23-mPDGRho)<cutRho) isRho=kTRUE;
 261   if(isRho){
 262     UInt_t pdg1[4]={321,211,211,211};
 263     Double_t mass1=InvMass(4,pdg1);
 264     if(TMath::Abs(mass1-mPDG)<cutMass) {isTrue=kTRUE;isGood=kTRUE;}
 265     if(isTrue) massD0bar[0]=mass1;
 266     isTrue=kFALSE;
 267     UInt_t pdg2[4]={211,321,211,211};
 268     Double_t mass2=InvMass(4,pdg2);
 269     if(TMath::Abs(mass2-mPDG)<cutMass) {isTrue=kTRUE;isGood=kTRUE;}
 270     if(isTrue) massD0[0]=mass2;
 271     isTrue=kFALSE;
 272   }
 273   return isGood;
 274 }