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 Bool_t AliHFEsignalCuts::IsCharmElectron(const TObject * const o) const {
 142   //
 143   // Check if mother is coming from Charm
 144   //
 145   if(!dynamic_cast<const AliVParticle *>(o)) return kFALSE;
 146   Int_t esources = GetElecSource(dynamic_cast<const AliVParticle *>(o));
 147   if(esources == AliHFEmcQA::kDirectCharm)  // 1: direct D->e
 148     return kTRUE;
 149   else
 150     return kFALSE;
 151 }
 152
 153 //____________________________________________________________
 154 Bool_t AliHFEsignalCuts::IsBeautyElectron(const TObject * const o) const {
 155   //
 156   // Check if mother is coming from Beauty
 157   //
 158   if(!dynamic_cast<const AliVParticle *>(o)) return kFALSE;
 159   Int_t esources = GetElecSource(dynamic_cast<const AliVParticle *>(o));
 160   if(esources == AliHFEmcQA::kDirectBeauty || esources == AliHFEmcQA::kBeautyCharm)  // 2: B->e 3: B->D->e
 161     return kTRUE;
 162   else
 163     return kFALSE;
 164 }
 165
 166 //____________________________________________________________
 167 Bool_t AliHFEsignalCuts::IsGammaElectron(const TObject * const o) const {
 168   //
 169   // Check for MC if the electron is coming from Gamma
 170   //
 171   if(!dynamic_cast<const AliVParticle *>(o)) return kFALSE;
 172   Int_t esources = GetElecSource(dynamic_cast<const AliVParticle *>(o));
 173   if(esources >= AliHFEmcQA::kGammaPi0 && esources <= AliHFEmcQA::kGammaRho0 )  // 4: conversion electrons
 174   //if(esources == AliHFEmcQA::kGammaPi0 || esources == AliHFEmcQA::kGammaEta || esources == AliHFEmcQA::kGammaOmega || esources == AliHFEmcQA::kGammaPhi || esources == AliHFEmcQA::kGammaEtaPrime || esources == AliHFEmcQA::kGammaRho0 )  // 4: conversion electrons
 175     return kTRUE;
 176   else
 177     return kFALSE;
 178 }
 179
 180 //____________________________________________________________
 181 Bool_t AliHFEsignalCuts::IsNonHFElectron(const TObject * const o) const {
 182   //
 183   // Check for MC if the electron is coming from NonHFE except for conversion
 184   //
 185   if(!dynamic_cast<const AliVParticle *>(o)) return kFALSE;
 186   Int_t esources = GetElecSource(dynamic_cast<const AliVParticle *>(o));
 187   if(esources == AliHFEmcQA:: kPi0 || esources == AliHFEmcQA::kEta || esources == AliHFEmcQA::kOmega || esources == AliHFEmcQA::kPhi || esources == AliHFEmcQA::kEtaPrime || esources == AliHFEmcQA::kRho0)  // 4: conversion electrons
 188     return kTRUE;
 189   else
 190     return kFALSE;
 191 }
 192
 193 //____________________________________________________________
 194 Bool_t AliHFEsignalCuts::IsJpsiElectron(const TObject * const o) const {
 195   //
 196   // Check if mother is coming from Charm
 197   //
 198   if(!dynamic_cast<const AliVParticle *>(o)) return kFALSE;
 199   Int_t esources = GetElecSource(dynamic_cast<const AliVParticle *>(o));
 200   if(esources == AliHFEmcQA::kJpsi)  // 5: J/psi->ee
 201     return kTRUE;
 202   else
 203     return kFALSE;
 204 }
 205
 206 //____________________________________________________________
 207 Bool_t AliHFEsignalCuts::IsB2JpsiElectron(const TObject * const o) const {
 208   //
 209   // Check if mother is coming from Charm
 210   //
 211   if(!dynamic_cast<const AliVParticle *>(o)) return kFALSE;
 212   Int_t esources = GetElecSource(dynamic_cast<const AliVParticle *>(o));
 213   if(esources == AliHFEmcQA::kB2Jpsi)  // 6: B->Jpsi->ee
 214     return kTRUE;
 215   else
 216     return kFALSE;
 217 }
 218
 219 //____________________________________________________________
 220 Bool_t AliHFEsignalCuts::IsKe3Electron(const TObject * const o) const {
 221   //
 222   // Check if mother is coming from Charm
 223   //
 224   if(!dynamic_cast<const AliVParticle *>(o)) return kFALSE;
 225   Int_t esources = GetElecSource(dynamic_cast<const AliVParticle *>(o));
 226   if(esources == AliHFEmcQA::kKe3)  // 7: K->e
 227     return kTRUE;
 228   else
 229     return kFALSE;
 230 }
 231
 232 /*
 233 //____________________________________________________________
 234 Bool_t AliHFEsignalCuts::IsCharmElectron(const TObject * const o) const {
 235   //
 236   // Check if mother is coming from Charm
 237   //
 238   if(TMath::Abs(GetTrackPDG(dynamic_cast<const AliVParticle *>(o))) != 11) return kFALSE;
 239   Int_t motherpdg = TMath::Abs(GetMotherPDG(dynamic_cast<const AliVParticle *>(o)));
 240   AliDebug(1, Form("Mother PDG %d\n", motherpdg));
 241
 242   if((motherpdg % 1000) / 100 == 4) return kTRUE;    // charmed meson, 3rd position in pdg code == 4
 243   if(motherpdg / 1000 == 4) return kTRUE;            // charmed baryon, 4th position in pdg code == 4
 244   AliDebug(1, "No Charm\n");
 245   return kFALSE;
 246 }
 247
 248 //____________________________________________________________
 249 Bool_t AliHFEsignalCuts::IsBeautyElectron(const TObject * const o) const {
 250   //
 251   // Check if mother is coming from Beauty
 252   //
 253   if(TMath::Abs(GetTrackPDG(dynamic_cast<const AliVParticle *>(o))) != 11) return kFALSE;
 254   Int_t motherpdg = TMath::Abs(GetMotherPDG(dynamic_cast<const AliVParticle *>(o)));
 255   AliDebug(1, Form("Mother PDG %d\n", motherpdg));
 256
 257   if((motherpdg % 1000) / 100 == 5) return kTRUE;   // beauty meson, 3rd position in pdg code == 5
 258   if(motherpdg / 1000 == 5) return kTRUE;           // beauty baryon, 4th position in pdg code == 5
 259   AliDebug(1, "No Beauty\n");
 260   return kFALSE;
 261 }
 262
 263 //____________________________________________________________
 264 Bool_t AliHFEsignalCuts::IsGammaElectron(const TObject * const o) const {
 265   //
 266   // Check for MC if the electron is coming from Gamma
 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!=22){
 273     AliDebug(1, "No Gamma");
 274     return kFALSE;
 275   } else {
 276     AliDebug(1, "Gamma");
 277     return kTRUE;
 278   }
 279 }
 280 */
 281
 282 //____________________________________________________________
 283 Int_t AliHFEsignalCuts::GetMotherPDG(const AliVParticle * const track) const {
 284   //
 285   // Get Mother Pdg code for reconstructed respectively MC tracks
 286   //
 287   TClass *type = track->IsA();
 288   if((!fMC && (type == AliESDtrack::Class())) || (!fAODArrayMCInfo && (type == AliAODTrack::Class()))){
 289     AliDebug(1, "No MC Event Available\n");
 290     return 0;
 291   }
 292   const AliVParticle *motherParticle = NULL, *mctrack = NULL;
 293   Int_t label = TMath::Abs(track->GetLabel());
 294   if(type == AliESDtrack::Class()){
 295     // Reconstructed track
 296     if(label) mctrack = fMC->GetTrack(TMath::Abs(label));
 297
 298   }
 299   else if(type == AliAODTrack::Class()) {
 300     // MCParticle
 301     if(label && label < fAODArrayMCInfo->GetEntriesFast())
 302       mctrack = (AliVParticle *) fAODArrayMCInfo->At(label);
 303   }
 304   else {
 305     mctrack=track;
 306   }
 307
 308   if(!mctrack) return 0;
 309
 310   Int_t motherPDG = 0;
 311   if(TString(mctrack->IsA()->GetName()).CompareTo("AliMCParticle") == 0){
 312     // case MC Particle
 313     const AliMCParticle *esdmctrack = dynamic_cast<const AliMCParticle *>(mctrack);
 314     if(esdmctrack) motherParticle = fMC->GetTrack(esdmctrack->Particle()->GetFirstMother());
 315     if(motherParticle){
 316       const AliMCParticle *esdmcmother = dynamic_cast<const AliMCParticle *>(motherParticle);
 317       if(esdmcmother) motherPDG = TMath::Abs(esdmcmother->Particle()->GetPdgCode());
 318     }
 319   } else {
 320     // case AODMCParticle
 321     const AliAODMCParticle *aodmctrack = dynamic_cast<const AliAODMCParticle *>(mctrack);
 322     if(aodmctrack) {
 323       if(aodmctrack->GetMother() && aodmctrack->GetMother() < fAODArrayMCInfo->GetEntriesFast())
 324       motherParticle =  (AliVParticle *) fAODArrayMCInfo->At(aodmctrack->GetMother());
 325     }
 326     if(motherParticle){
 327       const AliAODMCParticle *aodmcmother = dynamic_cast<const AliAODMCParticle *>(motherParticle);
 328       if(aodmcmother) motherPDG = TMath::Abs(aodmcmother->GetPdgCode());
 329     }
 330   }
 331   return motherPDG;
 332 }
 333
 334 //____________________________________________________________
 335 Int_t AliHFEsignalCuts::GetTrackPDG(const AliVParticle * const track) const {
 336   //
 337   // Return PDG code of a particle itself
 338   //
 339   TClass *type = track->IsA();
 340   if((!fMC && (type == AliESDtrack::Class())) || (!fAODArrayMCInfo && (type == AliAODTrack::Class()))){
 341     AliDebug(1, "No MC Event Available\n");
 342     return 0;
 343   }
 344   const AliVParticle *mctrack = NULL;
 345   Int_t label = TMath::Abs(track->GetLabel());
 346   if(type == AliESDtrack::Class()){
 347     // Reconstructed track
 348     if(label) mctrack = fMC->GetTrack(TMath::Abs(label));
 349
 350   }
 351   else if(type == AliAODTrack::Class()) {
 352     // MCParticle
 353     if(label && label < fAODArrayMCInfo->GetEntriesFast())
 354       mctrack =  (AliVParticle *) fAODArrayMCInfo->At(label);
 355   }
 356   else {
 357     mctrack=track;
 358   }
 359
 360   if(!mctrack) return 0;
 361
 362   TString mctype = mctrack->IsA()->GetName();
 363   Int_t trackPdg = 0;
 364   if(!mctype.CompareTo("AliMCParticle")){
 365     const AliMCParticle *esdmc = dynamic_cast<const AliMCParticle *>(mctrack);
 366     if(esdmc) trackPdg = esdmc->Particle()->GetPdgCode();
 367   } else {
 368     const AliAODMCParticle *aodmc = dynamic_cast< const AliAODMCParticle *>(mctrack);
 369     if(aodmc) trackPdg = aodmc->GetPdgCode();
 370   }
 371   return trackPdg;
 372 }
 373
 374 //____________________________________________________________
 375 Int_t AliHFEsignalCuts::GetElecSource(const AliVParticle * const track) const {
 376   //
 377   // Return PDG code of a particle itself
 378   //
 379
 380   TClass *type = track->IsA();
 381   if((!fMC && (type == AliESDtrack::Class())) || (!fAODArrayMCInfo && (type == AliAODTrack::Class()))){
 382     AliDebug(1, "No MC Event Available\n");
 383     return 0;
 384   }
 385   if(!fMCQA){
 386     AliDebug(1, "No MCQA Available\n");
 387     return 0;
 388   }
 389   if(!track){
 390     AliDebug(1, "Track not Available\n");
 391     return 0;
 392   }
 393
 394   const AliVParticle *mctrack = NULL;
 395   TParticle *mcpart = NULL;
 396   Int_t label = TMath::Abs(track->GetLabel());
 397   //AliMCParticle *esdmcmother = NULL;
 398   if(type == AliESDtrack::Class()){
 399     // Reconstructed track
 400     if(label) mctrack = fMC->GetTrack(TMath::Abs(label));
 401
 402   }
 403   else if(type == AliAODTrack::Class()) {
 404     // MCParticle
 405     if(label && label < fAODArrayMCInfo->GetEntriesFast())
 406       mctrack = (AliVParticle *) fAODArrayMCInfo->At(label);
 407   }
 408   else {
 409     mctrack=track;
 410   }
 411   if(!mctrack) return 0;
 412
 413   Int_t eSource = 0;
 414   if(mctrack->IsA() == AliMCParticle::Class()){
 415     const AliMCParticle *esdmc = dynamic_cast<const AliMCParticle *>(mctrack);
 416     if(esdmc){
 417       mcpart = esdmc->Particle();
 418       eSource=fMCQA->GetElecSource(mcpart);
 419 /* // considering secondary pions
 420       if(eSource>=AliHFEmcQA::kGammaPi0) {  // conversion electron, be careful with the enum odering
 421         Int_t glabel=TMath::Abs(esdmc->GetMother()); // gamma label
 422         if((esdmcmother= dynamic_cast<AliMCParticle *>(fMC->GetTrack(glabel)))){
 423           glabel=TMath::Abs(esdmcmother->GetMother()); // gamma's mother's label
 424           if((esdmcmother= dynamic_cast<AliMCParticle *>(fMC->GetTrack(glabel)))){
 425             if(glabel>fMC->GetNumberOfPrimaries()) eSource=AliHFEmcQA::kElse;
 426           }
 427         }
 428       }
 429       else if(eSource==AliHFEmcQA::kPi0 || (eSource>=AliHFEmcQA::kEta && eSource<=AliHFEmcQA::kRho0) ){ // nonHFE except for the conversion electron
 430         Int_t glabel=TMath::Abs(esdmc->GetMother());
 431         if((esdmcmother= dynamic_cast<AliMCParticle *>(fMC->GetTrack(glabel)))){
 432           if(glabel>fMC->GetNumberOfPrimaries()) eSource=AliHFEmcQA::kElse;
 433         }
 434       }
 435 */
 436     }
 437   } else {
 438     eSource=fMCQA->GetElecSource(mctrack);
 439   }
 440   return eSource;
 441 }