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