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