PWGHF/hfe/AliHFEsignalCuts.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, 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 // Signal cuts
  17 // Checks whether a particle (reconstructed or MC) is coming from MC Signal
  18 // For more information see implementation file
  19 //
  20 // Autor:
  21 //   Markus Fasel <M.Fasel@gsi.de>
  22 //
  23 #include <TClass.h>
  24 #include <TMath.h>
  25 #include <TParticle.h>
  26 #include <TString.h>
  27
  28 #include "AliAODTrack.h"
  29 #include "AliAODMCParticle.h"
  30 #include "AliESDtrack.h"
  31 #include "AliLog.h"
  32 #include "AliMCEvent.h"
  33 #include "AliMCParticle.h"
  34 #include "AliVParticle.h"
  35 #include "TClonesArray.h"
  36
  37 #include "AliHFEsignalCuts.h"
  38 #include "AliHFEmcQA.h"
  39
  40 ClassImp(AliHFEsignalCuts)
  41
  42 //____________________________________________________________
  43 AliHFEsignalCuts::AliHFEsignalCuts():
  44   AliAnalysisCuts(),
  45   fMC(NULL),
  46   fAODArrayMCInfo(NULL),
  47   fMCQA(NULL)
  48 {
  49   //
  50   // Dummy constructor
  51   //
  52 }
  53
  54 //____________________________________________________________
  55 AliHFEsignalCuts::AliHFEsignalCuts(const Char_t *name, const Char_t *title):
  56   AliAnalysisCuts(name, title),
  57   fMC(NULL),
  58   fAODArrayMCInfo(NULL),
  59   fMCQA(NULL)
  60 {
  61   //
  62   // Default constructor
  63   //
  64   fMCQA = new AliHFEmcQA;
  65   if(fMCQA) fMCQA->Init();
  66 }
  67
  68 //____________________________________________________________
  69 AliHFEsignalCuts::AliHFEsignalCuts(const AliHFEsignalCuts &ref):
  70   AliAnalysisCuts(ref),
  71   fMC(ref.fMC),
  72   fAODArrayMCInfo(NULL),
  73   fMCQA(ref.fMCQA)
  74 {
  75   //
  76   // Copy constructor
  77   //
  78 }
  79
  80 //____________________________________________________________
  81 AliHFEsignalCuts &AliHFEsignalCuts::operator=(const AliHFEsignalCuts &ref){
  82   //
  83   // Assignment operator
  84   //
  85   if(this != &ref){
  86     //fMC = ref.fMC;
  87     //fAODArrayMCInfo = ref.fAODArrayMCInfo;
  88     fMCQA = ref.fMCQA;
  89   }
  90   return *this;
  91 }
  92
  93 //____________________________________________________________
  94 AliHFEsignalCuts::~AliHFEsignalCuts(){
  95   //
  96   // Destructor
  97   //
  98   if(fMCQA) delete fMCQA;
  99 }
 100
 101 //____________________________________________________________
 102 void AliHFEsignalCuts::SetMCEvent(AliMCEvent *mc){
 103   //
 104   // Set mc event
 105   //
 106   fMC = mc;
 107   if(fMCQA) fMCQA->SetMCEvent(mc);
 108 }
 109
 110 //____________________________________________________________
 111 void AliHFEsignalCuts::SetMCAODInfo(TClonesArray *mcarray){
 112   //
 113   // Set mc array info
 114   //
 115   fAODArrayMCInfo = mcarray;
 116   if(fMCQA) fMCQA->SetMCArray(mcarray);
 117 }
 118
 119 //____________________________________________________________
 120 Bool_t AliHFEsignalCuts::IsSelected(TObject *o){
 121   //
 122   // Define signal as electron coming from charm or beauty
 123   // @TODO: Implement setter so that also either of them can be defined
 124   // as signal alone
 125
 126
 127   return IsCharmElectron(o) || IsBeautyElectron(o);
 128 /*
 129   //saving time?
 130   Int_t esources = GetElecSource(dynamic_cast<const AliVParticle *>(o));
 131   if(esources>0)printf("esources= %d\n",esources);
 132   if(esources == AliHFEmcQA::kDirectCharm || esources == AliHFEmcQA::kDirectBeauty || esources == AliHFEmcQA::kBeautyCharm)  // 1: direct D->e, 2: B->e 3: B->D->e
 133     return kTRUE;
 134   else
 135     return kFALSE;
 136 */
 137
 138 }
 139
 140 //____________________________________________________________
 141 AliHFEsignalCuts::ESignalSource_t AliHFEsignalCuts::GetSignalSource(const TObject *const o) const{
 142   //
 143   // Get source type of the electron
 144   //
 145   ESignalSource_t source = kOther;
 146   const AliVParticle *v = dynamic_cast<const AliVParticle *>(o);
 147   if(!v) return source;
 148   Int_t esources = GetElecSource(v);
 149   if(esources == AliHFEmcQA::kDirectCharm) source = kEleCharm;
 150   else if(esources == AliHFEmcQA::kDirectBeauty || esources == AliHFEmcQA::kBeautyCharm) source = kEleBeauty;
 151   else if(esources >= AliHFEmcQA::kGammaPi0 && esources <= AliHFEmcQA::kGammaRho0) source = kEleGamma;
 152   else if(esources == AliHFEmcQA:: kPi0 || esources == AliHFEmcQA::kEta || esources == AliHFEmcQA::kOmega || esources == AliHFEmcQA::kPhi || esources == AliHFEmcQA::kEtaPrime || esources == AliHFEmcQA::kRho0) source = kEleNonHFE;
 153   else if(esources == AliHFEmcQA::kJpsi) source = kEleJPsi;
 154   else if(esources == AliHFEmcQA::kB2Jpsi) source = kEleBtoJPsi;
 155   else if(esources == AliHFEmcQA::kKe3) source = kEleKe3;
 156   return source;
 157 }
 158
 159 /*****************************************
 160  *        Old legacy code                *
 161  *****************************************/
 162
 163 //____________________________________________________________
 164 Bool_t AliHFEsignalCuts::IsCharmElectronOld(const TObject * const o) const {
 165   //
 166   // Check if mother is coming from Charm
 167   //
 168   const AliVParticle *v = dynamic_cast<const AliVParticle *>(o);
 169   if(!v) return kFALSE;
 170   Int_t esources = GetElecSource(v);
 171   if(esources == AliHFEmcQA::kDirectCharm)  // 1: direct D->e
 172     return kTRUE;
 173   else
 174     return kFALSE;
 175 }
 176
 177 //____________________________________________________________
 178 Bool_t AliHFEsignalCuts::IsBeautyElectronOld(const TObject * const o) const {
 179   //
 180   // Check if mother is coming from Beauty
 181   //
 182   const AliVParticle *v = dynamic_cast<const AliVParticle *>(o);
 183   if(!v) return kFALSE;
 184   Int_t esources = GetElecSource(v);
 185   if(esources == AliHFEmcQA::kDirectBeauty || esources == AliHFEmcQA::kBeautyCharm)  // 2: B->e 3: B->D->e
 186     return kTRUE;
 187   else
 188     return kFALSE;
 189 }
 190
 191 //____________________________________________________________
 192 Bool_t AliHFEsignalCuts::IsGammaElectronOld(const TObject * const o) const {
 193   //
 194   // Check for MC if the electron is coming from Gamma
 195   //
 196   const AliVParticle *v = dynamic_cast<const AliVParticle *>(o);
 197   if(!v) return kFALSE;
 198   Int_t esources = GetElecSource(v);
 199   if(esources >= AliHFEmcQA::kGammaPi0 && esources <= AliHFEmcQA::kGammaRho0 )  // 4: conversion electrons
 200   //if(esources == AliHFEmcQA::kGammaPi0 || esources == AliHFEmcQA::kGammaEta || esources == AliHFEmcQA::kGammaOmega || esources == AliHFEmcQA::kGammaPhi || esources == AliHFEmcQA::kGammaEtaPrime || esources == AliHFEmcQA::kGammaRho0 )  // 4: conversion electrons
 201     return kTRUE;
 202   else
 203     return kFALSE;
 204 }
 205
 206 //____________________________________________________________
 207 Bool_t AliHFEsignalCuts::IsNonHFElectronOld(const TObject * const o) const {
 208   //
 209   // Check for MC if the electron is coming from NonHFE except for conversion
 210   //
 211   const AliVParticle *v = dynamic_cast<const AliVParticle *>(o);
 212   if(!v) return kFALSE;
 213   Int_t esources = GetElecSource(v);
 214   if(esources == AliHFEmcQA:: kPi0 || esources == AliHFEmcQA::kEta || esources == AliHFEmcQA::kOmega || esources == AliHFEmcQA::kPhi || esources == AliHFEmcQA::kEtaPrime || esources == AliHFEmcQA::kRho0)  // 4: conversion electrons
 215     return kTRUE;
 216   else
 217     return kFALSE;
 218 }
 219
 220 //____________________________________________________________
 221 Bool_t AliHFEsignalCuts::IsJpsiElectronOld(const TObject * const o) const {
 222   //
 223   // Check if mother is coming from Charm
 224   //
 225   const AliVParticle *v = dynamic_cast<const AliVParticle *>(o);
 226   if(!v) return kFALSE;
 227   Int_t esources = GetElecSource(v);
 228   if(esources == AliHFEmcQA::kJpsi)  // 5: J/psi->ee
 229     return kTRUE;
 230   else
 231     return kFALSE;
 232 }
 233
 234 //____________________________________________________________
 235 Bool_t AliHFEsignalCuts::IsB2JpsiElectronOld(const TObject * const o) const {
 236   //
 237   // Check if mother is coming from Charm
 238   //
 239   const AliVParticle *v = dynamic_cast<const AliVParticle *>(o);
 240   if(!v) return kFALSE;
 241   Int_t esources = GetElecSource(v);
 242   if(esources == AliHFEmcQA::kB2Jpsi)  // 6: B->Jpsi->ee
 243     return kTRUE;
 244   else
 245     return kFALSE;
 246 }
 247
 248 //____________________________________________________________
 249 Bool_t AliHFEsignalCuts::IsKe3ElectronOld(const TObject * const o) const {
 250   //
 251   // Check if mother is coming from Charm
 252   //
 253   const AliVParticle *v = dynamic_cast<const AliVParticle *>(o);
 254   if(!v) return kFALSE;
 255   Int_t esources = GetElecSource(v);
 256   if(esources == AliHFEmcQA::kKe3)  // 7: K->e
 257     return kTRUE;
 258   else
 259     return kFALSE;
 260 }
 261
 262 /*
 263 //____________________________________________________________
 264 Bool_t AliHFEsignalCuts::IsCharmElectron(const TObject * const o) const {
 265   //
 266   // Check if mother is coming from Charm
 267   //
 268   if(TMath::Abs(GetTrackPDG(dynamic_cast<const AliVParticle *>(o))) != 11) return kFALSE;
 269   Int_t motherpdg = TMath::Abs(GetMotherPDG(dynamic_cast<const AliVParticle *>(o)));
 270   AliDebug(1, Form("Mother PDG %d\n", motherpdg));
 271
 272   if((motherpdg % 1000) / 100 == 4) return kTRUE;    // charmed meson, 3rd position in pdg code == 4
 273   if(motherpdg / 1000 == 4) return kTRUE;            // charmed baryon, 4th position in pdg code == 4
 274   AliDebug(1, "No Charm\n");
 275   return kFALSE;
 276 }
 277
 278 //____________________________________________________________
 279 Bool_t AliHFEsignalCuts::IsBeautyElectron(const TObject * const o) const {
 280   //
 281   // Check if mother is coming from Beauty
 282   //
 283   if(TMath::Abs(GetTrackPDG(dynamic_cast<const AliVParticle *>(o))) != 11) return kFALSE;
 284   Int_t motherpdg = TMath::Abs(GetMotherPDG(dynamic_cast<const AliVParticle *>(o)));
 285   AliDebug(1, Form("Mother PDG %d\n", motherpdg));
 286
 287   if((motherpdg % 1000) / 100 == 5) return kTRUE;   // beauty meson, 3rd position in pdg code == 5
 288   if(motherpdg / 1000 == 5) return kTRUE;           // beauty baryon, 4th position in pdg code == 5
 289   AliDebug(1, "No Beauty\n");
 290   return kFALSE;
 291 }
 292
 293 //____________________________________________________________
 294 Bool_t AliHFEsignalCuts::IsGammaElectron(const TObject * const o) const {
 295   //
 296   // Check for MC if the electron is coming from Gamma
 297   //
 298   if(TMath::Abs(GetTrackPDG(dynamic_cast<const AliVParticle *>(o))) != 11) return kFALSE;
 299   Int_t motherpdg = TMath::Abs(GetMotherPDG(dynamic_cast<const AliVParticle *>(o)));
 300   AliDebug(1, Form("Mother PDG %d\n", motherpdg));
 301
 302   if(motherpdg!=22){
 303     AliDebug(1, "No Gamma");
 304     return kFALSE;
 305   } else {
 306     AliDebug(1, "Gamma");
 307     return kTRUE;
 308   }
 309 }
 310 */
 311
 312 //____________________________________________________________
 313 Int_t AliHFEsignalCuts::GetMotherPDG(const AliVParticle * const track) const {
 314   //
 315   // Get Mother Pdg code for reconstructed respectively MC tracks
 316   //
 317   TClass *type = track->IsA();
 318   //if((!fMC && (type == AliESDtrack::Class())) || (!fAODArrayMCInfo && (type == AliAODTrack::Class()))){
 319   //  AliDebug(1, "No MC Event Available\n");
 320   //  return 0;
 321   //}
 322   const AliVParticle *motherParticle = NULL, *mctrack = NULL;
 323   Int_t label = TMath::Abs(track->GetLabel());
 324   if(type == AliESDtrack::Class()){
 325     //
 326     if(!fMC) {
 327       AliDebug(1, "No MC Event Available\n");
 328       return 0;
 329     }
 330     // Reconstructed track
 331     if(label) mctrack = fMC->GetTrack(TMath::Abs(label));
 332
 333   }
 334   else if(type == AliAODTrack::Class()) {
 335     //
 336     if(!fAODArrayMCInfo) {
 337       AliDebug(1, "No MC Event Available\n");
 338       return 0;
 339     }
 340     // MCParticle
 341     if(label && label < fAODArrayMCInfo->GetEntriesFast())
 342       mctrack = (AliVParticle *) fAODArrayMCInfo->At(label);
 343   }
 344   else {
 345     mctrack=track;
 346   }
 347
 348   if(!mctrack) return 0;
 349
 350   Int_t motherPDG = 0;
 351   if(TString(mctrack->IsA()->GetName()).CompareTo("AliMCParticle") == 0){
 352     // case MC Particle
 353     const AliMCParticle *esdmctrack = dynamic_cast<const AliMCParticle *>(mctrack);
 354     if(esdmctrack) {
 355       if(!fMC) {
 356         AliDebug(1, "No MC Event Available\n");
 357         return 0;
 358       }
 359       motherParticle = fMC->GetTrack(esdmctrack->Particle()->GetFirstMother());
 360     }
 361     if(motherParticle){
 362       const AliMCParticle *esdmcmother = dynamic_cast<const AliMCParticle *>(motherParticle);
 363       if(esdmcmother) motherPDG = TMath::Abs(esdmcmother->Particle()->GetPdgCode());
 364     }
 365
 366   } else {
 367     // case AODMCParticle
 368     const AliAODMCParticle *aodmctrack = dynamic_cast<const AliAODMCParticle *>(mctrack);
 369     if(aodmctrack) {
 370       if(!fAODArrayMCInfo) {
 371         AliDebug(1, "No MC Event Available\n");
 372         return 0;
 373       }
 374       if(aodmctrack->GetMother() && aodmctrack->GetMother() < fAODArrayMCInfo->GetEntriesFast())
 375       motherParticle =  (AliVParticle *) fAODArrayMCInfo->At(aodmctrack->GetMother());
 376     }
 377     if(motherParticle){
 378       const AliAODMCParticle *aodmcmother = dynamic_cast<const AliAODMCParticle *>(motherParticle);
 379       if(aodmcmother) motherPDG = TMath::Abs(aodmcmother->GetPdgCode());
 380     }
 381   }
 382   return motherPDG;
 383 }
 384
 385 //____________________________________________________________
 386 Int_t AliHFEsignalCuts::GetTrackPDG(const AliVParticle * const track) const {
 387   //
 388   // Return PDG code of a particle itself
 389   //
 390
 391
 392   TClass *type = track->IsA();
 393   //if((!fMC && (type == AliESDtrack::Class())) || (!fAODArrayMCInfo && (type == AliAODTrack::Class()))){
 394   //  AliDebug(1, "No MC Event Available\n");
 395   //  return 0;
 396   //}
 397   const AliVParticle *mctrack = NULL;
 398   Int_t label = TMath::Abs(track->GetLabel());
 399   if(type == AliESDtrack::Class()){
 400     //
 401     if(!fMC) {
 402       AliDebug(1, "No MC Event Available\n");
 403       return 0;
 404     }
 405     // Reconstructed track
 406     if(label) mctrack = fMC->GetTrack(TMath::Abs(label));
 407
 408   }
 409   else if(type == AliAODTrack::Class()) {
 410     //
 411     if(!fAODArrayMCInfo) {
 412       AliDebug(1, "No MC Event Available\n");
 413       return 0;
 414     }
 415     // MCParticle
 416     if(label && label < fAODArrayMCInfo->GetEntriesFast())
 417       mctrack =  (AliVParticle *) fAODArrayMCInfo->At(label);
 418   }
 419   else {
 420     mctrack=track;
 421   }
 422
 423  if(!mctrack) return 0;
 424
 425   TString mctype = mctrack->IsA()->GetName();
 426   Int_t trackPdg = 0;
 427   if(!mctype.CompareTo("AliMCParticle")){
 428     const AliMCParticle *esdmc = dynamic_cast<const AliMCParticle *>(mctrack);
 429     if(esdmc) trackPdg = esdmc->Particle()->GetPdgCode();
 430   } else {
 431     const AliAODMCParticle *aodmc = dynamic_cast< const AliAODMCParticle *>(mctrack);
 432     if(aodmc) trackPdg = aodmc->GetPdgCode();
 433   }
 434   return trackPdg;
 435 }
 436
 437 //____________________________________________________________
 438 Int_t AliHFEsignalCuts::GetElecSource(const AliVParticle * const track) const {
 439   //
 440   // Return PDG code of a particle itself
 441   //
 442   if(!track){
 443     AliDebug(1, "Track not Available\n");
 444     return 0;
 445   }
 446
 447   TClass *type = track->IsA();
 448   //if((!fMC && (type == AliESDtrack::Class())) || (!fAODArrayMCInfo && (type == AliAODTrack::Class()))){
 449   //  AliDebug(1, "No MC Event Available\n");
 450   //  return 0;
 451   //}
 452   if(!fMCQA){
 453     AliDebug(1, "No MCQA Available\n");
 454     return 0;
 455   }
 456
 457
 458   const AliVParticle *mctrack = NULL;
 459   TParticle *mcpart = NULL;
 460   Int_t label = TMath::Abs(track->GetLabel());
 461   AliMCParticle *esdmcmother = NULL;
 462   if(type == AliESDtrack::Class()){
 463     //
 464     if(!fMC) {
 465       AliDebug(1, "No MC Event Available\n");
 466       return 0;
 467     }
 468     // Reconstructed track
 469     if(label) mctrack = fMC->GetTrack(TMath::Abs(label));
 470
 471   }
 472   else if(type == AliAODTrack::Class()) {
 473     //
 474     if(!fAODArrayMCInfo) {
 475       AliDebug(1, "No MC Event Available\n");
 476       return 0;
 477     }
 478     // MCParticle
 479     if(label && label < fAODArrayMCInfo->GetEntriesFast())
 480       mctrack = (AliVParticle *) fAODArrayMCInfo->At(label);
 481   }
 482   else {
 483     mctrack=track;
 484   }
 485   if(!mctrack) return 0;
 486
 487   Int_t eSource = 0;
 488   if(mctrack->IsA() == AliMCParticle::Class()){
 489     const AliMCParticle *esdmc = dynamic_cast<const AliMCParticle *>(mctrack);
 490     if(esdmc){
 491       mcpart = esdmc->Particle();
 492       eSource=fMCQA->GetElecSource(mcpart,kTRUE);
 493       // considering secondary pions
 494       if(type == AliESDtrack::Class()){
 495        if(eSource>=AliHFEmcQA::kGammaPi0 && eSource<=AliHFEmcQA::kGammaRho0) {  // conversion electron, be careful with the enum odering
 496         Int_t glabel=TMath::Abs(esdmc->GetMother()); // gamma label
 497         if((esdmcmother= dynamic_cast<AliMCParticle *>(fMC->GetTrack(glabel)))){
 498           glabel=TMath::Abs(esdmcmother->GetMother()); // gamma's mother's label
 499           if((esdmcmother= dynamic_cast<AliMCParticle *>(fMC->GetTrack(glabel)))){
 500             if(glabel>fMC->GetNumberOfPrimaries()) eSource=AliHFEmcQA::kScdryM;
 501           }
 502         }
 503        }
 504        else if(eSource==AliHFEmcQA::kPi0 || (eSource>=AliHFEmcQA::kEta && eSource<=AliHFEmcQA::kRho0) ){ // nonHFE except for the conversion electron
 505         Int_t glabel=TMath::Abs(esdmc->GetMother());
 506         if((esdmcmother= dynamic_cast<AliMCParticle *>(fMC->GetTrack(glabel)))){
 507           if(glabel>fMC->GetNumberOfPrimaries()) eSource=AliHFEmcQA::kScdryM;
 508         }
 509        }
 510       }
 511     }
 512   } else {
 513     eSource=fMCQA->GetElecSource(mctrack,kTRUE);
 514   }
 515   return eSource;
 516 }