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