PWGLF/SPECTRA/Kinks/AliResonanceKinkLikeSign.cxx

   1 /**************************************************************************
   2  * Author: Paraskevi Ganoti, University of Athens (pganoti@phys.uoa.gr)   *
   3  * Contributors are mentioned in the code where appropriate.              *
   4  *                                                                        *
   5  * Permission to use, copy, modify and distribute this software and its   *
   6  * documentation strictly for non-commercial purposes is hereby granted   *
   7  * without fee, provided that the above copyright notice appears in all   *
   8  * copies and that both the copyright notice and this permission notice   *
   9  * appear in the supporting documentation. The authors make no claims     *
  10  * about the suitability of this software for any purpose. It is          *
  11  * provided "as is" without express or implied warranty.                  *
  12  **************************************************************************/
  13
  14 //----------------------------------------------------------------------------------------------------------------
  15 //                        class AliResonanceKinkLikeSign
  16 //        Example of an analysis task for producing a like-sign background for resonances having at least one
  17 //        kaon-kink in their decay products.
  18 //        Background is calculated from a positive kaon kink and a negative track.
  19 //-----------------------------------------------------------------------------------------------------------------
  20
  21 #include "AliESDEvent.h"
  22 #include "TH1D.h"
  23 #include "TH2D.h"
  24 #include "TF1.h"
  25 #include "TDatabasePDG.h"
  26 #include "TLorentzVector.h"
  27 #include "AliESDkink.h"
  28
  29 #include "AliResonanceKinkLikeSign.h"
  30
  31 ClassImp(AliResonanceKinkLikeSign)
  32
  33 //________________________________________________________________________
  34 AliResonanceKinkLikeSign::AliResonanceKinkLikeSign(const char *name)
  35   : AliAnalysisTaskSE(name), fDebug(0), fListOfHistos(0), f1(0), f2(0), fPosKaonLikeSign(0), fLikeSignInvmassPt(0), fMaxNSigmaToVertex(0), fMinPtTrackCut(0), fMaxDCAxy(0), fMaxDCAzaxis(0),
  36 fMinTPCclusters(0),fMaxChi2PerTPCcluster(0), fMaxCov0(0), fMaxCov2(0), fMaxCov5(0) , fMaxCov9(0), fMaxCov14(0), fdaughter1pdg(0), fdaughter2pdg(0), fnbins(0), fnlowx(0), fnhighx(0), floweta(0), fuppereta(0), fminKinkRadius(0), fmaxKinkRadius(0), fminQt(0), fmaxQt(0), fptbins(0), flowpt(0), fupperpt(0), fmaxAbsEtaCut(0)
  37
  38 {
  39   // Constructor
  40
  41   // Define input and output slots here
  42   // Input slot #0 works with a TChain
  43   DefineOutput(1, TList::Class());
  44 }
  45
  46 //________________________________________________________________________
  47 void AliResonanceKinkLikeSign::UserCreateOutputObjects()
  48 {
  49   // Create histograms
  50   // Called once
  51
  52    f1=new TF1("f1","((atan([0]*[1]*(1.0/(sqrt((x^2)*(1.0-([1]^2))-([0]^2)*([1]^2))))))*180.)/[2]",1.1,10.0);
  53    f1->SetParameter(0,0.493677);
  54    f1->SetParameter(1,0.9127037);
  55    f1->SetParameter(2,TMath::Pi());
  56
  57    f2=new TF1("f2","((atan([0]*[1]*(1.0/(sqrt((x^2)*(1.0-([1]^2))-([0]^2)*([1]^2))))))*180.)/[2]",0.1,10.0);
  58    f2->SetParameter(0,0.13957018);
  59    f2->SetParameter(1,0.2731374);
  60    f2->SetParameter(2,TMath::Pi());
  61
  62    fPosKaonLikeSign=new TH1D("fPosKaonLikeSign"," ", fnbins, fnlowx, fnhighx);
  63    fLikeSignInvmassPt=new TH2D("fLikeSignInvmassPt"," ", fnbins, fnlowx, fnhighx, fptbins, flowpt, fupperpt);
  64
  65    fListOfHistos=new TList();
  66    fListOfHistos->Add(fPosKaonLikeSign);
  67    fListOfHistos->Add(fLikeSignInvmassPt);
  68
  69 }
  70
  71 //________________________________________________________________________
  72 void AliResonanceKinkLikeSign::UserExec(Option_t *)
  73 {
  74   // Main loop
  75   // Called for each event
  76
  77   AliVEvent *event = InputEvent();
  78   if (!event) {
  79      Printf("ERROR: Could not retrieve event");
  80      return;
  81   }
  82
  83   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
  84   if (!esd) {
  85      Printf("ERROR: Could not retrieve esd");
  86      return;
  87   }
  88
  89    Double_t daughter1Mass, daughter2Mass;
  90
  91    if (fdaughter1pdg==kKPlus)  {
  92      daughter1Mass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass();
  93      daughter2Mass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass();
  94    }
  95
  96    if (fdaughter1pdg!=kKPlus)  {
  97      daughter1Mass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass();
  98      daughter2Mass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass();
  99    }  // to ensure that daughter1Mass has always the kaon mass
 100
 101   if (!esd) {
 102     Printf("ERROR: esd not available");
 103     return;
 104   }
 105
 106   const AliESDVertex* vertex = GetEventVertex(esd);
 107   if(!vertex) return;
 108
 109   Double_t ptrackpos[3], ptrackneg[3];
 110
 111   TLorentzVector p4pos;
 112   TLorentzVector p4neg;
 113   TLorentzVector p4comb;
 114
 115   for (Int_t iTracks = 0; iTracks < esd->GetNumberOfTracks(); iTracks++) {
 116
 117     AliESDtrack* trackpos = esd->GetTrack(iTracks);
 118     if (!trackpos) {
 119       if (fDebug > 0) Printf("ERROR: Could not receive track %d", iTracks);
 120       continue;
 121     }
 122     if (trackpos->GetSign() < 0) continue;
 123
 124     AliExternalTrackParam *tpcTrackpos = (AliExternalTrackParam *)trackpos->GetTPCInnerParam();
 125     if(!tpcTrackpos) continue;
 126     ptrackpos[0]=tpcTrackpos->Px();
 127     ptrackpos[1]=tpcTrackpos->Py();
 128     ptrackpos[2]=tpcTrackpos->Pz();
 129
 130     Bool_t firstLevelAcceptPosTrack=IsAcceptedForKink(esd, vertex, trackpos);
 131     if(firstLevelAcceptPosTrack==kFALSE) continue;
 132
 133     TVector3 posTrackMom(ptrackpos[0],ptrackpos[1],ptrackpos[2]);
 134
 135     Int_t indexKinkPos=trackpos->GetKinkIndex(0);
 136     Bool_t posKaonKinkFlag=0;
 137     if(indexKinkPos<0) posKaonKinkFlag=IsKink(esd, indexKinkPos, posTrackMom);
 138     if(posKaonKinkFlag==0) continue;
 139     if(posKaonKinkFlag==1) p4pos.SetVectM(posTrackMom, daughter1Mass);
 140
 141     for (Int_t j=0; j<esd->GetNumberOfTracks(); j++) {
 142         if(iTracks==j) continue;
 143         AliESDtrack* trackneg=esd->GetTrack(j);
 144         if (!trackneg) {
 145           if (fDebug > 0) Printf("ERROR: Could not receive track %d", j);
 146           continue;
 147         }
 148         if (trackneg->GetSign() < 0) continue;
 149
 150         AliExternalTrackParam *tpcTrackneg = (AliExternalTrackParam *)trackneg->GetTPCInnerParam();
 151         if(!tpcTrackneg) continue;
 152         ptrackneg[0]=tpcTrackneg->Px();
 153         ptrackneg[1]=tpcTrackneg->Py();
 154         ptrackneg[2]=tpcTrackneg->Pz();
 155
 156         Bool_t firstLevelAcceptNegTrack=IsAcceptedForKink(esd, vertex, trackneg);
 157         if(firstLevelAcceptNegTrack==kFALSE) continue;
 158
 159         TVector3 negTrackMom(ptrackneg[0],ptrackneg[1],ptrackneg[2]);
 160
 161         Bool_t secondLevelAcceptNegTrack=IsAcceptedForTrack(esd, vertex, trackneg);
 162         if(secondLevelAcceptNegTrack==kFALSE) continue;
 163
 164         p4neg.SetVectM(negTrackMom, daughter2Mass);
 165
 166         p4comb=p4pos;
 167         p4comb+=p4neg;
 168
 169         if(p4comb.Vect().Pt()<=fMinPtTrackCut) continue;
 170
 171         if((TMath::Abs(p4pos.Vect().Eta())<fmaxAbsEtaCut)&&(TMath::Abs(p4neg.Vect().Eta())<fmaxAbsEtaCut)&&(p4comb.Vect().Eta()<fmaxAbsEtaCut)) {
 172
 173           fPosKaonLikeSign->Fill(p4comb.M());
 174           fLikeSignInvmassPt->Fill(p4comb.M(), p4comb.Vect().Pt());
 175
 176         }
 177
 178       } //inner track loop
 179
 180   } //outer track loop
 181
 182     PostData(1, fListOfHistos);
 183 }
 184
 185 //________________________________________________________________________
 186 void AliResonanceKinkLikeSign::Terminate(Option_t *)
 187 {
 188
 189 }
 190
 191 //____________________________________________________________________//
 192
 193 Float_t AliResonanceKinkLikeSign::GetSigmaToVertex(AliESDtrack* esdTrack) const {
 194   // Calculates the number of sigma to the vertex.
 195
 196   Float_t b[2];
 197   Float_t bRes[2];
 198   Float_t bCov[3];
 199
 200     esdTrack->GetImpactParameters(b,bCov);
 201
 202   if (bCov[0]<=0 || bCov[2]<=0) {
 203     //AliDebug(1, "Estimated b resolution lower or equal zero!");
 204     bCov[0]=0; bCov[2]=0;
 205   }
 206   bRes[0] = TMath::Sqrt(bCov[0]);
 207   bRes[1] = TMath::Sqrt(bCov[2]);
 208
 209   if (bRes[0] == 0 || bRes[1] ==0) return -1;
 210
 211   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 212
 213   if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
 214
 215   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 216
 217   return d;
 218 }
 219
 220 //____________________________________________________________________//
 221
 222 const AliESDVertex* AliResonanceKinkLikeSign::GetEventVertex(const AliESDEvent* esd) const
 223
 224 {
 225   // Get the vertex
 226
 227   const AliESDVertex* vertex = esd->GetPrimaryVertex();
 228
 229   if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
 230   else
 231   {
 232      vertex = esd->GetPrimaryVertexSPD();
 233      if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
 234      else
 235      return 0;
 236   }
 237 }
 238
 239 //________________________________________________________________________
 240
 241  Bool_t AliResonanceKinkLikeSign::IsAcceptedForKink(AliESDEvent *localesd,
 242             const AliESDVertex *localvertex, AliESDtrack* localtrack) {
 243    // Apply the selections for each kink
 244
 245   Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
 246   Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};  //The impact parameters and their covariance.
 247   Double_t dca3D = 0.0;
 248
 249   AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)localtrack->GetTPCInnerParam();
 250   if(!tpcTrack) {
 251     gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
 252     dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
 253     cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
 254   }
 255   else {
 256     gPt = tpcTrack->Pt();
 257     gPx = tpcTrack->Px();
 258     gPy = tpcTrack->Py();
 259     gPz = tpcTrack->Pz();
 260     tpcTrack->PropagateToDCA(localvertex,
 261                localesd->GetMagneticField(),100.,dca,cov);
 262   }
 263
 264   if(GetSigmaToVertex(localtrack) > fMaxNSigmaToVertex) {
 265       if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a %lf sigmas to vertex TPC (max. requested: %lf)",   GetSigmaToVertex(localtrack),fMaxNSigmaToVertex);
 266       return kFALSE;
 267   }
 268
 269   if(TMath::Abs(dca[0]) > fMaxDCAxy) {
 270       if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a value of dca(xy) (TPC) of %lf (max. requested: %lf)", TMath::Abs(dca[0]), fMaxDCAxy);
 271       return kFALSE;
 272   }
 273
 274   if(TMath::Abs(dca[1]) > fMaxDCAzaxis) {
 275       if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a value of dca(z) of %lf (max. requested: %lf)", TMath::Abs(dca[1]), fMaxDCAzaxis);
 276       return kFALSE;
 277   }
 278
 279   if(gPt <= fMinPtTrackCut) {
 280       if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a min value of pt of %lf (min. requested: %lf)", gPt, fMinPtTrackCut);
 281       return kFALSE;
 282   }
 283
 284   return kTRUE;
 285 }
 286
 287 //________________________________________________________________________
 288 Bool_t AliResonanceKinkLikeSign::IsAcceptedForTrack(AliESDEvent *localesd,                                                                                                                                          const AliESDVertex *localvertex, AliESDtrack *localtrack) {
 289    // Apply the selections for each track
 290
 291   Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
 292   Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};  //The impact parameters and their covariance.
 293   Double_t dca3D = 0.0;
 294
 295   AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)localtrack->GetTPCInnerParam();
 296   if(!tpcTrack) {
 297     gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
 298     dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
 299     cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
 300   }
 301   else {
 302     gPt = tpcTrack->Pt();
 303     gPx = tpcTrack->Px();
 304     gPy = tpcTrack->Py();
 305     gPz = tpcTrack->Pz();
 306     tpcTrack->PropagateToDCA(localvertex,
 307                localesd->GetMagneticField(),100.,dca,cov);
 308   }
 309
 310   Int_t fcls[200];
 311   Int_t nClustersTPC=localtrack->GetTPCclusters(fcls);
 312   Float_t chi2perTPCcluster=-1.0;
 313   if(nClustersTPC!=0) chi2perTPCcluster=(localtrack->GetTPCchi2())/Float_t(nClustersTPC);
 314
 315   Double_t extCov[15];
 316   localtrack->GetExternalCovariance(extCov);
 317
 318   if((localtrack->GetStatus() & AliESDtrack::kTPCrefit) == 0) {
 319       if (fDebug > 1) Printf("IsAccepted: Track rejected because of no refited in TPC");
 320       return kFALSE;
 321   }
 322
 323   if(nClustersTPC < fMinTPCclusters) {
 324       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of nclusters (TPC) of %d (min. requested: %d)", nClustersTPC, fMinTPCclusters);
 325       return kFALSE;
 326   }
 327
 328   if(chi2perTPCcluster > fMaxChi2PerTPCcluster) {
 329       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of chi2perTPCcluster of %lf (max. requested: %lf)", chi2perTPCcluster, fMaxChi2PerTPCcluster);
 330       return kFALSE;
 331   }
 332
 333   if(extCov[0] > fMaxCov0) {
 334       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[0] of %lf (max. requested: %lf)", extCov[0], fMaxCov0);
 335       return kFALSE;
 336   }
 337
 338   if(extCov[2] > fMaxCov2) {
 339       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[2] of %lf (max. requested: %lf)", extCov[2], fMaxCov2);
 340       return kFALSE;
 341   }
 342
 343   if(extCov[5] > fMaxCov5) {
 344       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[5] of %lf (max. requested: %lf)", extCov[5], fMaxCov5);
 345       return kFALSE;
 346   }
 347
 348   if(extCov[9] > fMaxCov9) {
 349       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[9] of %lf (max. requested: %lf)", extCov[9], fMaxCov9);
 350       return kFALSE;
 351   }
 352
 353   if(extCov[14] > fMaxCov14) {
 354       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[14] of %lf (max. requested: %lf)", extCov[14], fMaxCov14);
 355       return kFALSE;
 356   }
 357
 358   return kTRUE;
 359
 360 }
 361
 362 //________________________________________________________________________
 363 Bool_t AliResonanceKinkLikeSign::IsKink(AliESDEvent *localesd, Int_t kinkIndex, TVector3 trackMom)
 364 {
 365    // Test some kinematical criteria for each kink
 366
 367          AliESDkink *kink=localesd->GetKink(TMath::Abs(kinkIndex)-1);
 368          const TVector3 motherMfromKink(kink->GetMotherP());
 369          const TVector3 daughterMKink(kink->GetDaughterP());
 370          Float_t qt=kink->GetQt();
 371
 372          Double_t maxDecAngKmu=f1->Eval(motherMfromKink.Mag(),0.,0.,0.);
 373          Double_t maxDecAngpimu=f2->Eval(motherMfromKink.Mag(),0.,0.,0.);
 374
 375          Float_t kinkAngle=TMath::RadToDeg()*kink->GetAngle(2);
 376
 377          Float_t energyDaughterMu=TMath::Sqrt(daughterMKink.Mag()*daughterMKink.Mag()+0.105658*0.105658);
 378          Float_t p1XM= motherMfromKink.Px();
 379          Float_t p1YM= motherMfromKink.Py();
 380          Float_t p1ZM= motherMfromKink.Pz();
 381          Float_t p2XM= daughterMKink.Px();
 382          Float_t p2YM= daughterMKink.Py();
 383          Float_t p2ZM= daughterMKink.Pz();
 384          Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM)));
 385          Double_t invariantMassKmu= TMath::Sqrt((energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-motherMfromKink.Mag()*motherMfromKink.Mag());
 386
 387          if((kinkAngle>maxDecAngpimu)&&(qt>fminQt)&&(qt<fmaxQt)&&((kink->GetR()>fminKinkRadius)&&(kink->GetR()<fmaxKinkRadius))&&(TMath::Abs(trackMom.Eta())<fmaxAbsEtaCut)&&(invariantMassKmu<0.6)) {
 388
 389            if(trackMom.Mag()<=1.1) {
 390                 return kTRUE;
 391            }
 392            else
 393            if (kinkAngle<maxDecAngKmu) {
 394                 return kTRUE;
 395            }
 396          }
 397          return kFALSE;
 398 }