PWGHF/correlationHF/AliDxHFEParticleSelectionMCEl.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   AliDxHFEParticleSelectionMCEl.cxx
  20 /// @author Hege Erdal, Matthias Richter
  21 /// @date   2012-07-19
  22 /// @brief  MC El selection for D0-HFE correlation
  23 ///
  24
  25 #include "AliDxHFEParticleSelectionMCEl.h"
  26 #include "AliVParticle.h"
  27 #include "AliLog.h"
  28 #include "THnSparse.h"
  29 #include "AliAODMCParticle.h"
  30 #include "TH1F.h"
  31 #include "TAxis.h"
  32 #include "AliAODTrack.h"
  33 #include "AliReducedParticle.h"
  34 #include <iostream>
  35 #include <cerrno>
  36 #include <memory>
  37
  38 using namespace std;
  39 using std::vector;
  40
  41
  42 /// ROOT macro for the implementation of ROOT specific class methods
  43 ClassImp(AliDxHFEParticleSelectionMCEl)
  44
  45 AliDxHFEParticleSelectionMCEl::AliDxHFEParticleSelectionMCEl(const char* opt)
  46   : AliDxHFEParticleSelectionEl(opt)
  47   , fMCTools()
  48   , fPDGnotMCElectron(NULL)
  49   , fOriginMother(0)
  50   , fResultMC(0)
  51   , fUseKine(kFALSE)
  52   , fMotherPDGs(0)
  53   , fUseMCReco(kFALSE)
  54   , fSelectionStep(AliDxHFEParticleSelectionEl::kNoCuts)
  55 {
  56   // constructor
  57   //
  58   //
  59   //
  60   //
  61
  62   fMCTools.~AliDxHFEToolsMC();
  63
  64   ParseArguments(opt);
  65
  66   // This is also checked in AddTasks, but just for safety!
  67   if(fUseMCReco && fUseKine) AliFatal("CAN'T SET BOTH usekine AND elmcreco AT THE SAME TIME");
  68
  69   // TODO: argument scan, build tool options accordingly
  70   // e.g. set mc mode first/last, skip control histograms
  71   TString toolopt("pdg=11 mc-last");
  72   if(fUseKine) toolopt+=" usekine";
  73   new (&fMCTools) AliDxHFEToolsMC(toolopt);
  74 }
  75
  76 //at the moment not used, keep for now (to be used when plotting)
  77 const char* AliDxHFEParticleSelectionMCEl::fgkPDGMotherBinLabels[]={
  78   "d",
  79   "u",
  80   "s",
  81   "c",
  82   "b",
  83   "gluon",
  84   "gamma",
  85   "#pi^{0}",
  86   "#eta",
  87   "proton",
  88   "others"
  89 };
  90
  91
  92 const char*  AliDxHFEParticleSelectionMCEl::fgkPDGBinLabels[]={
  93   "positron",
  94   "electron",
  95   "#mu+",
  96   "#mu-",
  97   "#pi+",
  98   "#pi-",
  99   "K+",
 100   "K-",
 101   "proton",
 102   "antiproton",
 103   "others"
 104 };
 105
 106 AliDxHFEParticleSelectionMCEl::~AliDxHFEParticleSelectionMCEl()
 107 {
 108   // destructor
 109
 110   if(fPDGnotMCElectron){
 111     delete fPDGnotMCElectron;
 112     fPDGnotMCElectron=NULL;
 113   }
 114
 115 }
 116
 117 int AliDxHFEParticleSelectionMCEl::Init()
 118 {
 119   //
 120   // init function
 121   //
 122
 123   // Particles considered HFE background. can be expanded
 124   fMotherPDGs.push_back(AliDxHFEToolsMC::kPDGpi0);
 125   fMotherPDGs.push_back(AliDxHFEToolsMC::kPDGeta);
 126   fMotherPDGs.push_back(AliDxHFEToolsMC::kPDGgamma);
 127   fMotherPDGs.push_back(AliDxHFEToolsMC::kPDGJpsi);
 128
 129   int iResult=0;
 130   iResult=AliDxHFEParticleSelectionEl::Init();
 131   if (iResult<0) return iResult;
 132
 133   // Histo containing PDG of track which was not MC truth electron
 134   // TODO: Add to list in baseclass
 135   fPDGnotMCElectron=CreateControlHistogram("fPDGnotMCElectron","PDG of track not MC truth electron",AliDxHFEToolsMC::kNofPDGLabels,fgkPDGBinLabels);
 136   AddControlObject(fPDGnotMCElectron);
 137   return 0;
 138 }
 139
 140 THnSparse* AliDxHFEParticleSelectionMCEl::DefineTHnSparse()
 141 {
 142   //
 143   // Defines the THnSparse.
 144
 145   const int thnSize = 4;
 146   InitTHnSparseArray(thnSize);
 147   const double Pi=TMath::Pi();
 148   TString name;
 149   name.Form("%s info", GetName());
 150
 151   //                           0    1      2     3
 152   //                           Pt   Phi   Eta   mother
 153   int    thnBins[thnSize] = { 1000,  200, 500,    15  };
 154   double thnMin [thnSize] = {    0,    0, -1.,  -1.5  };
 155   double thnMax [thnSize] = {  100, 2*Pi,  1.,  13.5  };
 156   const char* thnNames[thnSize]={
 157     "Pt",
 158     "Phi",
 159     "Eta",
 160     "Mother", //bin==-1: Not MC truth electron
 161   };
 162
 163   return CreateControlTHnSparse(name,thnSize,thnBins,thnMin,thnMax,thnNames);
 164 }
 165
 166 int AliDxHFEParticleSelectionMCEl::FillParticleProperties(AliVParticle* p, Double_t* data, int dimension) const
 167 {
 168   // fill the data array from the particle data
 169   if (!data) return -EINVAL;
 170   AliAODTrack *track=(AliAODTrack*)p;
 171   if (!track) return -ENODATA;
 172   int i=0;
 173   if (dimension!=GetDimTHnSparse()) {
 174     // TODO: think about filling only the available data and throwing a warning
 175     return -ENOSPC;
 176   }
 177   data[i++]=track->Pt();
 178   data[i++]=track->Phi();
 179   data[i++]=track->Eta();
 180   data[i++]=fOriginMother;
 181
 182   return i;
 183 }
 184
 185 int AliDxHFEParticleSelectionMCEl::IsSelected(AliVParticle* p, const AliVEvent* pEvent)
 186 {
 187   /// overloaded from AliDxHFEParticleSelection: check particle
 188   /// H: Have changed function. Now doing particle selection first, then run MC over
 189   /// selected tracks. Could configure it to be configurable, but not sure if it
 190   /// is needed.
 191   /// Result from normal track selection is returned, result from MC is stored in
 192   /// THnSparse.
 193
 194   int iResult=0;
 195   if (!p || !pEvent){
 196     return -EINVAL;
 197   }
 198   fOriginMother=-1;
 199
 200   if(!fUseKine && !fUseMCReco){
 201   // step 1:
 202   // optional MC selection before the particle selection
 203   if (fMCTools.MCFirst() && (iResult=CheckMC(p, pEvent))==0) {
 204     // histograming?
 205     return iResult;
 206   }
 207
 208   // step 2 or 1, depending on sequence:
 209   // normal particle selection
 210   iResult=AliDxHFEParticleSelectionEl::IsSelected(p, pEvent);
 211   if (fMCTools.MCFirst() || iResult==0) return iResult;
 212
 213
 214   }
 215
 216   // step 2, only executed if MC check is last
 217   // optional MC selection after the particle selection
 218   // result stored to be filled into THnSparse
 219   // TODO: strictly speaken the particles should be rejected
 220   // if not mc selected, however skip this for the moment, because of
 221   // the logic outside
 222   // This line will run always when running directly over the stack or over MC reconstructed tracks
 223   // For MC reconstructed tracks, should consider adding more constrictions on the tracks,
 224   // maybe do the track selection first (which means could call AliDxHFEParticleSelectionEl::IsSelected()
 225   // and don't do PID
 226
 227   if(fUseMCReco && ( fSelectionStep > AliDxHFEParticleSelectionEl::kNoCuts )){
 228     iResult=AliDxHFEParticleSelectionEl::IsSelected(p, pEvent);
 229     if(iResult ==0) return iResult;
 230   }
 231   fResultMC=CheckMC(p, pEvent);
 232
 233   if(fUseKine || fUseMCReco)
 234     return fResultMC;
 235
 236   return iResult;
 237
 238 }
 239
 240 int AliDxHFEParticleSelectionMCEl::CheckMC(AliVParticle* p, const AliVEvent* pEvent)
 241 {
 242   /// check if MC criteria are fulfilled
 243   if (!p || !pEvent){
 244     return -EINVAL;
 245   }
 246   int iResult=0;
 247
 248   if (!fMCTools.IsInitialized() && (iResult=fMCTools.InitMCParticles(pEvent))<0) {
 249     // TODO: message? but has to be filtered in order to avoid message flood
 250     return 0; // no meaningful filtering on mc possible
 251   }
 252
 253   if(fUseKine){
 254     // If run on kinematical level, check if particle is electron (only ones who are relevant)
 255     // and if pass IsPhysicalPrimary()
 256     Int_t test = fMCTools.CheckMCParticle(p);
 257     if(test==1) return 0;
 258     // Add additional contraints on the electrons?
 259     // remove delta e? also remove E<300MeV?
 260     // Should also mark dalitz decay and gamma conversion..
 261     AliAODMCParticle *mcp=dynamic_cast<AliAODMCParticle*>(p);
 262     if(!mcp->IsPhysicalPrimary()) return 0;
 263
 264   }
 265
 266   int pdgParticle=-1;
 267   if (!fUseKine && fMCTools.RejectByPDG(p,false, &pdgParticle)) {
 268     // rejected by pdg
 269     // TODO: Move this to fMCTools???? Can this be part of the statistics in the MC class?
 270     fPDGnotMCElectron->Fill(fMCTools.MapPDGLabel(pdgParticle));
 271     return 0;
 272   }
 273
 274   int pdgMother=0;
 275   // Find PDG of first mother
 276   pdgMother=fMCTools.FindMotherPDG(p,AliDxHFEToolsMC::kGetFirstMother);
 277
 278   // Check if first mother is counted as background
 279   Bool_t isNotBackground=fMCTools.RejectByPDG(pdgMother,fMotherPDGs);
 280
 281   if(!isNotBackground){
 282     // Set fOriginMother if mother is counted as background
 283     // TODO: Could this be done in a more elegant way?
 284     switch(pdgMother){
 285     case(AliDxHFEToolsMC::kPDGpi0): fOriginMother=AliDxHFEToolsMC::kNrOrginMother; break;
 286     case(AliDxHFEToolsMC::kPDGeta): fOriginMother=AliDxHFEToolsMC::kNrOrginMother+1; break;
 287     case(AliDxHFEToolsMC::kPDGgamma): fOriginMother=AliDxHFEToolsMC::kNrOrginMother+2;break;
 288     case(AliDxHFEToolsMC::kPDGJpsi): fOriginMother=AliDxHFEToolsMC::kNrOrginMother+3;break;
 289     }
 290   }
 291   else{
 292     // If loop over Stack, also checks if First mother is a HF meson
 293     Bool_t isHFmeson =fMCTools.TestMotherHFMeson(TMath::Abs(pdgMother));
 294
 295     if(isHFmeson){
 296       // If first mother is a HF meson, loops back to find
 297       // original quark + if there was a gluon. Result is
 298       // stored in fOriginMother
 299       pdgMother=fMCTools.FindMotherPDG(p,AliDxHFEToolsMC::kGetOriginMother);
 300       fOriginMother=fMCTools.GetOriginMother();
 301     }
 302     else{
 303       fOriginMother=AliDxHFEToolsMC::kNrOrginMother+4;
 304     }
 305
 306   }
 307
 308   return 1;
 309 }
 310
 311 void AliDxHFEParticleSelectionMCEl::Clear(const char* option)
 312 {
 313   /// clear internal memory
 314   fMCTools.Clear(option);
 315 }
 316
 317 AliVParticle *AliDxHFEParticleSelectionMCEl::CreateParticle(AliVParticle* track)
 318 {
 319
 320   AliReducedParticle *part = new AliReducedParticle(track->Eta(), track->Phi(), track->Pt(),track->Charge(),fOriginMother);
 321
 322   return part;
 323
 324 }
 325
 326 int AliDxHFEParticleSelectionMCEl::ParseArguments(const char* arguments)
 327 {
 328   // parse arguments and set internal flags
 329   TString strArguments(arguments);
 330   auto_ptr<TObjArray> tokens(strArguments.Tokenize(" "));
 331   if (!tokens.get()) return 0;
 332
 333   AliInfo(strArguments);
 334   TIter next(tokens.get());
 335   TObject* token;
 336   while ((token=next())) {
 337     TString argument=token->GetName();
 338     if (argument.BeginsWith("usekine") ){
 339       fUseKine=kTRUE;
 340       continue;
 341     }
 342     if (argument.BeginsWith("elmcreco")){
 343       fUseMCReco=kTRUE;
 344       if(argument.BeginsWith("elmcreco=")){
 345         argument.ReplaceAll("elmcreco=", "");
 346         if(argument.CompareTo("aftertrackcuts")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kHFEcutsTPC;
 347         if(argument.CompareTo("aftertofpid")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kPIDTOF;
 348         if(argument.CompareTo("afterfullpid")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kPIDTOFTPC;
 349
 350         AliDxHFEParticleSelectionEl::SetFinalCutStep(fSelectionStep);
 351
 352       }
 353       continue;
 354     }
 355     // forwarding of single argument works, unless key-option pairs separated
 356     // by blanks are introduced
 357     AliDxHFEParticleSelection::ParseArguments(argument);
 358   }
 359
 360   return 0;
 361 }