PWGHF/correlationHF/AliDxHFEToolsMC.cxx

   1 // $Id$
   2
   3 //**************************************************************************
   4 //* This file is property of and copyright by the ALICE Project            *
   5 //* ALICE Experiment at CERN, All rights reserved.                         *
   6 //*                                                                        *
   7 //* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *
   8 //*                  Hege Erdal       <hege.erdal@gmail.com>               *
   9 //*                                                                        *
  10 //* Permission to use, copy, modify and distribute this software and its   *
  11 //* documentation strictly for non-commercial purposes is hereby granted   *
  12 //* without fee, provided that the above copyright notice appears in all   *
  13 //* copies and that both the copyright notice and this permission notice   *
  14 //* appear in the supporting documentation. The authors make no claims     *
  15 //* about the suitability of this software for any purpose. It is          *
  16 //* provided "as is" without express or implied warranty.                  *
  17 //**************************************************************************
  18
  19 /// @file   AliDxHFEToolsMC.cxx
  20 /// @author Hege Erdal, Matthias Richter
  21 /// @date   2012-07-19
  22 /// @brief  Common Tools for MC particle selection
  23 ///
  24
  25 #include "AliDxHFEToolsMC.h"
  26 #include "AliAODMCParticle.h"
  27 #include "AliVEvent.h"
  28 #include "AliVParticle.h"
  29 #include "AliLog.h"
  30 #include "TObjArray.h"
  31 #include "TString.h"
  32 #include "TH1D.h"
  33 #include <iostream>
  34 #include <cerrno>
  35 #include <memory>
  36
  37 using namespace std;
  38
  39 /// ROOT macro for the implementation of ROOT specific class methods
  40 ClassImp(AliDxHFEToolsMC)
  41
  42 AliDxHFEToolsMC::AliDxHFEToolsMC(const char* option)
  43   : fSequence(kMCLast)
  44   , fMCParticles(NULL)
  45   , fPDGs()
  46   , fMotherPDGs()
  47   , fHistPDG(NULL)
  48   , fHistPDGMother(NULL)
  49   , fOriginMother(kOriginNone)
  50   , fMClabel(-1)
  51   , fNrMCParticles(-1)
  52   , fUseKine(kFALSE)
  53 {
  54   // constructor
  55   //
  56   //
  57   //
  58   //
  59   Init(option);
  60 }
  61
  62 const char*  AliDxHFEToolsMC::fgkPDGBinLabels[]={
  63   "positron",
  64   "electron",
  65   "#mu+",
  66   "#mu-",
  67   "#pi+",
  68   "#pi-",
  69   "K+",
  70   "K-",
  71   "proton",
  72   "antiproton",
  73   "others"
  74 };
  75
  76 const char*  AliDxHFEToolsMC::fgkPDGMotherBinLabels[]={
  77   "d",
  78   "u",
  79   "s",
  80   "c",
  81   "b",
  82   "gluon",
  83   "gamma",
  84   "#pi^{0}",
  85   "#eta",
  86   "proton",
  87   "others"
  88 };
  89
  90 const char*  AliDxHFEToolsMC::fgkStatisticsBinLabels[]={
  91   "all",   // all MC particles
  92   "found", // selected particles with correct MC
  93   "fake"   // selected particles without corresponding MC
  94 };
  95
  96 AliDxHFEToolsMC::~AliDxHFEToolsMC()
  97 {
  98   // destructor
  99   if (fHistPDG) delete fHistPDG;
 100   fHistPDG=NULL;
 101   if (fHistPDGMother) delete fHistPDGMother;
 102   fHistPDGMother=NULL;
 103 }
 104
 105 int AliDxHFEToolsMC::Init(const char* option)
 106 {
 107   // initialize according to options
 108   TString strOption(option);
 109   bool bControlHist=true;
 110   std::auto_ptr<TObjArray> tokens(strOption.Tokenize(" "));
 111   if (tokens.get() && tokens->GetEntriesFast()>0) {
 112     for (int itoken=0; itoken<tokens->GetEntriesFast(); itoken++) {
 113       if (tokens->At(itoken)==NULL) continue;
 114       TString arg=tokens->At(itoken)->GetName();
 115       const char* key="";
 116       key="pdg=";
 117       if (arg.BeginsWith(key)) {
 118         arg.ReplaceAll(key, "");
 119         fPDGs.push_back(arg.Atoi());
 120       }
 121       key="mother-pdg=";
 122       if (arg.BeginsWith(key)) {
 123         arg.ReplaceAll(key, "");
 124         fMotherPDGs.push_back(arg.Atoi());
 125       }
 126       key="control-hist=";
 127       if (arg.BeginsWith(key)) {
 128         arg.ReplaceAll(key, "");
 129         bControlHist=arg.CompareTo("off")!=0;
 130       }
 131       key="mc-first";
 132       if (arg.BeginsWith(key)) {
 133         fSequence=kMCFirst;
 134       }
 135       key="mc-last";
 136       if (arg.BeginsWith(key)) {
 137         fSequence=kMCLast;
 138       }
 139       key="usekine";
 140       if(arg.BeginsWith(key)) {
 141         printf("AliDxHFEToolsMC::Init()  Using Kinematical\n");
 142         fUseKine=kTRUE;
 143       }
 144     }
 145   }
 146
 147   if (bControlHist) {
 148     fHistPDG=CreateControlHistogram("histPDG",
 149                                     "pdg code of selected particle",
 150                                     sizeof(fgkPDGBinLabels)/sizeof(const char*),
 151                                     fgkPDGBinLabels);
 152     fHistPDGMother=CreateControlHistogram("histPDGMother",
 153                                           "pdg code of first mother of selected particle",
 154                                           sizeof(fgkPDGMotherBinLabels)/sizeof(fgkPDGMotherBinLabels[0]),
 155                                           fgkPDGMotherBinLabels);
 156   }
 157   return 0;
 158 }
 159
 160 int AliDxHFEToolsMC::InitMCParticles(const AliVEvent* pEvent)
 161 {
 162
 163   // init MC info from event object
 164   if (!pEvent) return -EINVAL;
 165
 166   // TODO: choose branch name depending on VEvent type; configurable?
 167   TString branchname(AliAODMCParticle::StdBranchName());
 168   TObject* o=pEvent->FindListObject(branchname);
 169   if (!o) {
 170     AliWarningClass(Form("can not find MC info '%s' in event of type '%s'", branchname.Data(), pEvent->ClassName()));
 171     return -ENOENT;
 172   }
 173   fMCParticles = dynamic_cast<TObjArray*>(o);
 174   if (!fMCParticles) {
 175     AliWarningClass(Form("ignoring MC info '%s' of wrong type '%s', expecting TObjArray", branchname.Data(), o->ClassName()));
 176     return -ENODATA;
 177   }
 178   fNrMCParticles=fMCParticles->GetEntriesFast();
 179   return 0;
 180 }
 181
 182 bool AliDxHFEToolsMC::RejectByPDG(int pdg, const vector<int> &list) const
 183 {
 184   // check if pdg should be rejected, particle is not rejected
 185   // if it is in the list, returns always false if list is empty
 186   if (list.size()==0) return false;
 187   for (vector<int>::const_iterator i=list.begin();
 188        i!=list.end(); i++) {
 189     if (*i==pdg) return false;
 190   }
 191   return true;
 192 }
 193
 194 bool AliDxHFEToolsMC::RejectByPDG(AliVParticle* p, bool doStatistics, int* pdgParticleResult)
 195 {
 196   // check if pdg should be rejected
 197   // always false if not pdg list is initialized
 198   if (!p) return false;
 199   Int_t MClabel = p->GetLabel();
 200   if(MClabel<0) {
 201     return true;
 202   }
 203   int pdgPart=-1;
 204   // TODO: there might be more types of particles to be checked
 205   AliAODMCParticle* aodmcp=0;
 206   aodmcp=dynamic_cast<AliAODMCParticle*>(fMCParticles->At(MClabel));
 207   if (aodmcp)
 208     pdgPart=TMath::Abs(aodmcp->GetPdgCode());
 209   if (pdgPart<0) return 0;
 210
 211   if (pdgParticleResult)
 212     *pdgParticleResult=pdgPart;
 213
 214   bool bReject=RejectByPDG(pdgPart, fPDGs);
 215   if (doStatistics && fHistPDG) {
 216     // TODO: think about histogramming mode, e.g. histogramming of rejected particles?
 217     if (!bReject) {
 218       fHistPDG->Fill(MapPDGLabel(p->PdgCode()));
 219     }
 220   }
 221   return bReject;
 222 }
 223
 224 bool AliDxHFEToolsMC::RejectByMotherPDG(AliVParticle* p, bool doStatistics)
 225 {
 226   // check if pdg should be rejected by mother
 227   // always false if not mother pdg list is initialized
 228   // H: think maybe this is particle specific, and should be moved to PartSelMCEl
 229   if (!p) return false;
 230   int motherpdg=FindPdgOriginMother(p);
 231   // TODO: This should be tweaked. Want to
 232   bool bReject=RejectByPDG(motherpdg, fMotherPDGs);
 233   if (doStatistics && fHistPDGMother) {
 234     // TODO: think about histogramming mode, e.g. histogramming of rejected particles?
 235     if (!bReject) {
 236       fHistPDGMother->Fill(MapPDGMotherLabel(p->PdgCode()));
 237     }
 238   }
 239   return bReject;
 240 }
 241
 242 int AliDxHFEToolsMC::FindMotherPDG(AliVParticle* p, bool bReturnFirstMother)
 243 {
 244   // Will find and return pdg of the mother, either first or loop down to the
 245   // initial quark
 246
 247   // To reset fOriginMother.
 248   fOriginMother=kOriginNone;
 249
 250   if (!p) return false;
 251   int motherpdg=FindPdgOriginMother(p, bReturnFirstMother);
 252
 253   return motherpdg;
 254 }
 255
 256 TH1* AliDxHFEToolsMC::CreateControlHistogram(const char* name,
 257                                              const char* title,
 258                                              int nBins,
 259                                              const char** binLabels) const
 260 {
 261   /// create control histogram
 262   std::auto_ptr<TH1> h(new TH1D(name, title, nBins, -0.5, nBins-0.5));
 263   if (!h.get()) return NULL;
 264   for (int iLabel=0; iLabel<nBins; iLabel++) {
 265     h->GetXaxis()->SetBinLabel(iLabel+1, binLabels[iLabel]);
 266   }
 267
 268   return h.release();
 269 }
 270
 271 int AliDxHFEToolsMC::MapPDGLabel(int pdg) const
 272 {
 273   /// mapping of pdg code to enum
 274   switch (pdg) {
 275   case kPDGelectron  : return kPDGLabelElectron;
 276   case -kPDGelectron : return kPDGLabelPositron;
 277   case kPDGmuon      : return kPDGLabelMuPlus;
 278   case -kPDGmuon     : return kPDGLabelMuMinus;
 279   case kPDGpion      : return kPDGLabelPiPlus;
 280   case -kPDGpion     : return kPDGLabelPiMinus;
 281   case kPDGkaon      : return kPDGLabelKPlus;
 282   case -kPDGkaon     : return kPDGLabelKMinus;
 283   case kPDGproton    : return kPDGLabelProton;
 284   case -kPDGproton   : return kPDGLabelAntiproton;
 285   default:
 286     return kPDGLabelOthers;
 287   }
 288 }
 289
 290 int AliDxHFEToolsMC::MapPDGMotherLabel(int pdg) const
 291 {
 292   /// mapping of pdg code to enum
 293   switch (pdg) {
 294   case kPDGd     : return kPDGMotherLabelD;
 295   case kPDGu     : return kPDGMotherLabelU;
 296   case kPDGs     : return kPDGMotherLabelS;
 297   case kPDGc     : return kPDGMotherLabelC;
 298   case kPDGb     : return kPDGMotherLabelB;
 299   case kPDGgluon : return kPDGMotherLabelGluon;
 300   case kPDGgamma : return kPDGMotherLabelGamma;
 301   case kPDGpi0   : return kPDGMotherLabelPi0;
 302   case kPDGeta   : return kPDGMotherLabelEta;
 303   case kPDGproton: return kPDGMotherLabelProton;
 304   default:
 305     return kPDGLabelOthers;
 306   }
 307 }
 308
 309 int AliDxHFEToolsMC::FindPdgOriginMother(AliVParticle* p, bool bReturnFirstMother)
 310 {
 311   // Return the pgd of original mother particle
 312   // TODO: need also to have specific for D0, electron etc
 313   // for instance to mark when you have gluon, charm or beauty
 314   // among the mothers. Or maybe this will be the same anyway?
 315   // TODO: implement tests on origin, if charm/beauty quark and if
 316   // they came from gluon. use booleans to set this which can be accessed from
 317   // outside? Something like fSequence.
 318
 319   if (!p) return kPDGnone;
 320
 321   Int_t imother=-1;
 322   Int_t MClabel=0;
 323
 324   // Either using MClabel set from outside (needed for Dmesons), or find it
 325   if(fMClabel<0){
 326     MClabel = p->GetLabel();
 327     if(MClabel<0){
 328       return kPDGnone;
 329     }
 330   }
 331   else MClabel=fMClabel;
 332
 333   // try different classes, unfortunately there is no common base class
 334   AliAODMCParticle* aodmcp=0;
 335   if(fUseKine)
 336     aodmcp=dynamic_cast<AliAODMCParticle*>(p);
 337   else
 338     aodmcp=dynamic_cast<AliAODMCParticle*>(fMCParticles->At(MClabel));
 339
 340   if (!aodmcp) {
 341     return kPDGnone;
 342   }
 343   imother = aodmcp->GetMother();
 344
 345   // Absolute or +/- on pgd ???
 346   Int_t pdg=TMath::Abs(aodmcp->GetPdgCode());
 347
 348   if (imother<0){
 349     // also check this particle
 350     CheckOriginMother(pdg);
 351     return pdg;
 352   }
 353
 354   if (!fMCParticles->At(imother)) {
 355     AliErrorClass(Form("no mc particle with label %d", imother));
 356     return pdg;
 357   }
 358
 359   AliAODMCParticle * mother=dynamic_cast<AliAODMCParticle*>(fMCParticles->At(imother));
 360   if (!mother) {
 361     AliErrorClass(Form("mc mother particle of wrong class type %s", fMCParticles->At(imother)->ClassName()));
 362     return pdg;
 363   }
 364
 365   if(bReturnFirstMother){
 366     return mother->GetPdgCode();
 367   }
 368
 369   CheckOriginMother(pdg);
 370
 371   if (mother->GetPdgCode()==kPDGproton){
 372     // why? - H: This is the proton, can't get further back
 373     // To be discussed whether to do this a different way
 374     return pdg;
 375   }
 376
 377   //Reset fMClabel to find pdg of mother if looping on D
 378   fMClabel=-1;
 379   pdg=FindPdgOriginMother(mother);
 380
 381   return pdg;
 382 }
 383
 384 void AliDxHFEToolsMC::CheckOriginMother(int pdg)
 385 {
 386
 387   // Checking if the particle is a quark or gluon and setting fOriginMother accordingly.
 388   // Right now only check on quark. Need also check on whether it comes from D or B meson?
 389
 390   switch(pdg){
 391   case(kPDGc):
 392     fOriginMother = kOriginCharm; break;
 393   case(kPDGb):
 394     fOriginMother = kOriginBeauty; break;
 395   case(kPDGgluon):
 396     if(fOriginMother==kOriginCharm) fOriginMother=kOriginGluonCharm;
 397     else if(fOriginMother==kOriginBeauty) fOriginMother=kOriginGluonBeauty;
 398     else fOriginMother=kOriginGluon;
 399     break;
 400   case(kPDGd):
 401     if(!TestIfHFquark(fOriginMother))
 402       fOriginMother=kOriginDown;
 403     break;
 404   case(kPDGu):
 405     if(!TestIfHFquark(fOriginMother))
 406       fOriginMother=kOriginUp;
 407     break;
 408   case(kPDGs):
 409     if(!TestIfHFquark(fOriginMother))
 410       fOriginMother=kOriginStrange;
 411     break;
 412   }
 413 }
 414
 415 Bool_t AliDxHFEToolsMC::TestIfHFquark(int origin)
 416 {
 417
 418   // Checking if particle has been marked as charm/beauty quark
 419
 420   Bool_t test=kFALSE;
 421   switch(origin){
 422   case(kOriginCharm):
 423     test=kTRUE; break;
 424   case(kOriginBeauty):
 425     test=kTRUE; break;
 426   case(kOriginGluonCharm):
 427     test=kTRUE; break;
 428   case(kOriginGluonBeauty):
 429     test=kTRUE; break;
 430   }
 431   return test;
 432 }
 433
 434 Bool_t AliDxHFEToolsMC::TestMotherHFMeson(int pdg)
 435 {
 436
 437   // Checking if the pdg corresponds to a HF meson (D or B meson)
 438
 439   Bool_t isD = kFALSE;
 440   Bool_t isB = kFALSE;
 441   if((pdg>=400 && pdg <500) || (pdg>=4000 && pdg<5000 ))
 442     isD = kTRUE;
 443   if((pdg>=500 && pdg <600) || (pdg>=5000 && pdg<6000 ))
 444     {isD = kFALSE; isB = kTRUE;}
 445
 446   if(isD || isB) return kTRUE;
 447   else return kFALSE;
 448
 449 }
 450
 451
 452
 453
 454 void AliDxHFEToolsMC::Clear(const char* /*option*/)
 455 {
 456   // clear internal memory
 457   fMCParticles=NULL;
 458   fNrMCParticles=-1;
 459 }
 460
 461
 462 int AliDxHFEToolsMC::CheckMCParticle(AliVParticle* p){
 463
 464   // Checks if MC particle is desired particle
 465
 466   AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(p);
 467   if (!mcPart) {
 468     AliInfoClass("MC Particle not found in tree, skipping");
 469     return -1;
 470   }
 471
 472   Int_t PDG =TMath::Abs(mcPart->PdgCode());
 473   int bReject=RejectByPDG(PDG, fPDGs);
 474
 475   return bReject;
 476
 477 }