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 {
  53   // constructor
  54   //
  55   //
  56   //
  57   //
  58   Init(option);
  59 }
  60
  61 const char*  AliDxHFEToolsMC::fgkPDGBinLabels[]={
  62   "positron",
  63   "electron",
  64   "#mu+",
  65   "#mu-",
  66   "#pi+",
  67   "#pi-",
  68   "K+",
  69   "K-",
  70   "proton",
  71   "antiproton",
  72   "others"
  73 };
  74
  75 const char*  AliDxHFEToolsMC::fgkPDGMotherBinLabels[]={
  76   "d",
  77   "u",
  78   "s",
  79   "c",
  80   "b",
  81   "gluon",
  82   "gamma",
  83   "#pi^{0}",
  84   "#eta",
  85   "proton",
  86   "others"
  87 };
  88
  89 const char*  AliDxHFEToolsMC::fgkStatisticsBinLabels[]={
  90   "all",   // all MC particles
  91   "found", // selected particles with correct MC
  92   "fake"   // selected particles without corresponding MC
  93 };
  94
  95 AliDxHFEToolsMC::~AliDxHFEToolsMC()
  96 {
  97   // destructor
  98   if (fHistPDG) delete fHistPDG;
  99   fHistPDG=NULL;
 100   if (fHistPDGMother) delete fHistPDGMother;
 101   fHistPDGMother=NULL;
 102 }
 103
 104 int AliDxHFEToolsMC::Init(const char* option)
 105 {
 106   // initialize according to options
 107   TString strOption(option);
 108   bool bControlHist=true;
 109   std::auto_ptr<TObjArray> tokens(strOption.Tokenize(" "));
 110   if (tokens.get() && tokens->GetEntriesFast()>0) {
 111     for (int itoken=0; itoken<tokens->GetEntriesFast(); itoken++) {
 112       if (tokens->At(itoken)==NULL) continue;
 113       TString arg=tokens->At(itoken)->GetName();
 114       const char* key="";
 115       key="pdg=";
 116       if (arg.BeginsWith(key)) {
 117         arg.ReplaceAll(key, "");
 118         fPDGs.push_back(arg.Atoi());
 119       }
 120       key="mother-pdg=";
 121       if (arg.BeginsWith(key)) {
 122         arg.ReplaceAll(key, "");
 123         fMotherPDGs.push_back(arg.Atoi());
 124       }
 125       key="control-hist=";
 126       if (arg.BeginsWith(key)) {
 127         arg.ReplaceAll(key, "");
 128         bControlHist=arg.CompareTo("off")!=0;
 129       }
 130       key="mc-first";
 131       if (arg.BeginsWith(key)) {
 132         fSequence=kMCFirst;
 133       }
 134       key="mc-last";
 135       if (arg.BeginsWith(key)) {
 136         fSequence=kMCLast;
 137       }
 138     }
 139   }
 140
 141   if (bControlHist) {
 142     fHistPDG=CreateControlHistogram("histPDG",
 143                                     "pdg code of selected particle",
 144                                     sizeof(fgkPDGBinLabels)/sizeof(const char*),
 145                                     fgkPDGBinLabels);
 146     fHistPDGMother=CreateControlHistogram("histPDGMother",
 147                                           "pdg code of first mother of selected particle",
 148                                           sizeof(fgkPDGMotherBinLabels)/sizeof(fgkPDGMotherBinLabels[0]),
 149                                           fgkPDGMotherBinLabels);
 150   }
 151   return 0;
 152 }
 153
 154 int AliDxHFEToolsMC::InitMCParticles(const AliVEvent* pEvent)
 155 {
 156
 157   // init MC info from event object
 158   if (!pEvent) return -EINVAL;
 159
 160   // TODO: choose branch name depending on VEvent type; configurable?
 161   TString branchname(AliAODMCParticle::StdBranchName());
 162   TObject* o=pEvent->FindListObject(branchname);
 163   if (!o) {
 164     AliWarningClass(Form("can not find MC info '%s' in event of type '%s'", branchname.Data(), pEvent->ClassName()));
 165     return -ENOENT;
 166   }
 167   fMCParticles = dynamic_cast<TObjArray*>(o);
 168   if (!fMCParticles) {
 169     AliWarningClass(Form("ignoring MC info '%s' of wrong type '%s', expecting TObjArray", branchname.Data(), o->ClassName()));
 170     return -ENODATA;
 171   }
 172   fNrMCParticles=fMCParticles->GetEntriesFast();
 173   return 0;
 174 }
 175
 176 bool AliDxHFEToolsMC::RejectByPDG(int pdg, const vector<int> &list) const
 177 {
 178   // check if pdg should be rejected, particle is not rejected
 179   // if it is in the list, returns always false if list is empty
 180   if (list.size()==0) return false;
 181   for (vector<int>::const_iterator i=list.begin();
 182        i!=list.end(); i++) {
 183     if (*i==pdg) return false;
 184   }
 185   return true;
 186 }
 187
 188 bool AliDxHFEToolsMC::RejectByPDG(AliVParticle* p, bool doStatistics, int* pdgParticleResult)
 189 {
 190   // check if pdg should be rejected
 191   // always false if not pdg list is initialized
 192   if (!p) return false;
 193   Int_t MClabel = p->GetLabel();
 194   if(MClabel<0) return 0;
 195
 196   int pdgPart=-1;
 197   // TODO: there might be more types of particles to be checked
 198   AliAODMCParticle* aodmcp=0;
 199   aodmcp=dynamic_cast<AliAODMCParticle*>(fMCParticles->At(MClabel));
 200   if (aodmcp)
 201     pdgPart=TMath::Abs(aodmcp->GetPdgCode());
 202   if (pdgPart<0) return 0;
 203
 204   if (pdgParticleResult)
 205     *pdgParticleResult=pdgPart;
 206
 207   bool bReject=RejectByPDG(pdgPart, fPDGs);
 208   if (doStatistics && fHistPDG) {
 209     // TODO: think about histogramming mode, e.g. histogramming of rejected particles?
 210     if (!bReject) {
 211       fHistPDG->Fill(MapPDGLabel(p->PdgCode()));
 212     }
 213   }
 214   return bReject;
 215 }
 216
 217 bool AliDxHFEToolsMC::RejectByMotherPDG(AliVParticle* p, bool doStatistics)
 218 {
 219   // check if pdg should be rejected by mother
 220   // always false if not mother pdg list is initialized
 221   // H: think maybe this is particle specific, and should be moved to PartSelMCEl
 222   if (!p) return false;
 223   int motherpdg=FindPdgOriginMother(p);
 224   // TODO: This should be tweaked. Want to
 225   bool bReject=RejectByPDG(motherpdg, fMotherPDGs);
 226   if (doStatistics && fHistPDGMother) {
 227     // TODO: think about histogramming mode, e.g. histogramming of rejected particles?
 228     if (!bReject) {
 229       fHistPDGMother->Fill(MapPDGMotherLabel(p->PdgCode()));
 230     }
 231   }
 232   return bReject;
 233 }
 234
 235 int AliDxHFEToolsMC::FindMotherPDG(AliVParticle* p, bool bReturnFirstMother)
 236 {
 237   // Will find and return pdg of the mother, either first or loop down to the
 238   // initial quark
 239
 240   // To reset fOriginMother.
 241   fOriginMother=kOriginNone;
 242
 243   if (!p) return false;
 244   int motherpdg=FindPdgOriginMother(p, bReturnFirstMother);
 245
 246   return motherpdg;
 247 }
 248
 249
 250 TH1* AliDxHFEToolsMC::CreateControlHistogram(const char* name,
 251                                              const char* title,
 252                                              int nBins,
 253                                              const char** binLabels) const
 254 {
 255   /// create control histogram
 256   std::auto_ptr<TH1> h(new TH1D(name, title, nBins, -0.5, nBins-0.5));
 257   if (!h.get()) return NULL;
 258   for (int iLabel=0; iLabel<nBins; iLabel++) {
 259     h->GetXaxis()->SetBinLabel(iLabel+1, binLabels[iLabel]);
 260   }
 261
 262   return h.release();
 263 }
 264
 265 int AliDxHFEToolsMC::MapPDGLabel(int pdg) const
 266 {
 267   /// mapping of pdg code to enum
 268   switch (pdg) {
 269   case kPDGelectron  : return kPDGLabelElectron;
 270   case -kPDGelectron : return kPDGLabelPositron;
 271   case kPDGmuon      : return kPDGLabelMuPlus;
 272   case -kPDGmuon     : return kPDGLabelMuMinus;
 273   case kPDGpion      : return kPDGLabelPiPlus;
 274   case -kPDGpion     : return kPDGLabelPiMinus;
 275   case kPDGkaon      : return kPDGLabelKPlus;
 276   case -kPDGkaon     : return kPDGLabelKMinus;
 277   case kPDGproton    : return kPDGLabelProton;
 278   case -kPDGproton   : return kPDGLabelAntiproton;
 279   default:
 280     return kPDGLabelOthers;
 281   }
 282 }
 283
 284 int AliDxHFEToolsMC::MapPDGMotherLabel(int pdg) const
 285 {
 286   /// mapping of pdg code to enum
 287   switch (pdg) {
 288   case kPDGd     : return kPDGMotherLabelD;
 289   case kPDGu     : return kPDGMotherLabelU;
 290   case kPDGs     : return kPDGMotherLabelS;
 291   case kPDGc     : return kPDGMotherLabelC;
 292   case kPDGb     : return kPDGMotherLabelB;
 293   case kPDGgluon : return kPDGMotherLabelGluon;
 294   case kPDGgamma : return kPDGMotherLabelGamma;
 295   case kPDGpi0   : return kPDGMotherLabelPi0;
 296   case kPDGeta   : return kPDGMotherLabelEta;
 297   case kPDGproton: return kPDGMotherLabelProton;
 298   default:
 299     return kPDGLabelOthers;
 300   }
 301 }
 302
 303 int AliDxHFEToolsMC::FindPdgOriginMother(AliVParticle* p, bool bReturnFirstMother)
 304 {
 305   // Return the pgd of original mother particle
 306   // TODO: need also to have specific for D0, electron etc
 307   // for instance to mark when you have gluon, charm or beauty
 308   // among the mothers. Or maybe this will be the same anyway?
 309   // TODO: implement tests on origin, if charm/beauty quark and if
 310   // they came from gluon. use booleans to set this which can be accessed from
 311   // outside? Something like fSequence.
 312
 313   if (!p) return kPDGnone;
 314
 315   Int_t imother=-1;
 316   Int_t MClabel=0;
 317
 318   // Either using MClabel set from outside (needed for Dmesons), or find it
 319   if(fMClabel<0){
 320     MClabel = p->GetLabel();
 321     if(MClabel<0){
 322       return kPDGnone;
 323     }
 324   }
 325   else MClabel=fMClabel;
 326
 327   // try different classes, unfortunately there is no common base class
 328   AliAODMCParticle* aodmcp=0;
 329   aodmcp=dynamic_cast<AliAODMCParticle*>(fMCParticles->At(MClabel));
 330   if (!aodmcp) {
 331     return kPDGnone;
 332   }
 333   imother = aodmcp->GetMother();
 334
 335   // Absolute or +/- on pgd ???
 336   Int_t pdg=TMath::Abs(aodmcp->GetPdgCode());
 337
 338   if (imother<0){
 339     // also check this particle
 340     CheckOriginMother(pdg);
 341     return pdg;
 342   }
 343
 344   if (!fMCParticles->At(imother)) {
 345     AliErrorClass(Form("no mc particle with label %d", imother));
 346     return pdg;
 347   }
 348
 349   AliAODMCParticle * mother=dynamic_cast<AliAODMCParticle*>(fMCParticles->At(imother));
 350   if (!mother) {
 351     AliErrorClass(Form("mc mother particle of wrong class type %s", fMCParticles->At(imother)->ClassName()));
 352     return pdg;
 353   }
 354
 355   if(bReturnFirstMother){
 356     return mother->GetPdgCode();
 357   }
 358
 359   CheckOriginMother(pdg);
 360
 361   if (mother->GetPdgCode()==kPDGproton){
 362     // why? - H: This is the proton, can't get further back
 363     // To be discussed whether to do this a different way
 364     return pdg;
 365   }
 366
 367   //Reset fMClabel to find pdg of mother if looping on D
 368   fMClabel=-1;
 369   pdg=FindPdgOriginMother(mother);
 370
 371   return pdg;
 372 }
 373
 374 void AliDxHFEToolsMC::CheckOriginMother(int pdg)
 375 {
 376
 377   // Checking if the particle is a quark or gluon and setting fOriginMother accordingly.
 378   // Right now only check on quark. Need also check on whether it comes from D or B meson?
 379
 380   switch(pdg){
 381   case(kPDGc):
 382     fOriginMother = kOriginCharm; break;
 383   case(kPDGb):
 384     fOriginMother = kOriginBeauty; break;
 385   case(kPDGgluon):
 386     if(fOriginMother==kOriginCharm) fOriginMother=kOriginGluonCharm;
 387     else if(fOriginMother==kOriginBeauty) fOriginMother=kOriginGluonBeauty;
 388     else fOriginMother=kOriginGluon;
 389     break;
 390   case(kPDGd):
 391     if(!TestIfHFquark(fOriginMother))
 392       fOriginMother=kOriginDown;
 393     break;
 394   case(kPDGu):
 395     if(!TestIfHFquark(fOriginMother))
 396       fOriginMother=kOriginUp;
 397     break;
 398   case(kPDGs):
 399     if(!TestIfHFquark(fOriginMother))
 400       fOriginMother=kOriginStrange;
 401     break;
 402   }
 403 }
 404
 405 Bool_t AliDxHFEToolsMC::TestIfHFquark(int origin)
 406 {
 407
 408   // Checking if particle has been marked as charm/beauty quark
 409
 410   Bool_t test=kFALSE;
 411   switch(origin){
 412   case(kOriginCharm):
 413     test=kTRUE; break;
 414   case(kOriginBeauty):
 415     test=kTRUE; break;
 416   case(kOriginGluonCharm):
 417     test=kTRUE; break;
 418   case(kOriginGluonBeauty):
 419     test=kTRUE; break;
 420   }
 421   return test;
 422 }
 423
 424 void AliDxHFEToolsMC::Clear(const char* /*option*/)
 425 {
 426   // clear internal memory
 427   fMCParticles=NULL;
 428   fNrMCParticles=-1;
 429 }
 430