PWGCF/EBYE/PIDFluctuation/task/AliEbyEPidRatioEffCont.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: ALICE Offline.                                                 *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 //=========================================================================//
  18 //             AliEbyE Analysis for Particle Ratio Fluctuation             //
  19 //                   Deepika Rathee  | Satyajit Jena                       //
  20 //                   drathee@cern.ch | sjena@cern.ch                       //
  21 //                  Date: Wed Jul  9 18:38:30 CEST 2014                    //
  22 //          New approch to find particle ratio to reduce memory            //
  23 //                             (Test Only)                                 //
  24 //=========================================================================//
  25
  26 #include "TMath.h"
  27 #include "TAxis.h"
  28
  29 #include "AliESDEvent.h"
  30 #include "AliStack.h"
  31 #include "AliMCEvent.h"
  32 #include "AliESDtrackCuts.h"
  33 #include "AliAODEvent.h"
  34 #include "AliAODMCParticle.h"
  35
  36 #include "AliEbyEPidRatioEffCont.h"
  37
  38 using namespace std;
  39
  40
  41 ClassImp(AliEbyEPidRatioEffCont)
  42
  43 //________________________________________________________________________
  44 AliEbyEPidRatioEffCont::AliEbyEPidRatioEffCont() :
  45   AliEbyEPidRatioBase("EffCont", "EffCont"),
  46   fLabelsRec(NULL),
  47   fHnEffMc(NULL),
  48   fHnContMc(NULL),
  49   fHnEffRec(NULL),
  50   fHnContRec(NULL) {
  51   // Constructor
  52
  53   AliLog::SetClassDebugLevel("AliEbyEPidRatioEffCont",10);
  54 }
  55
  56 //________________________________________________________________________
  57 AliEbyEPidRatioEffCont::~AliEbyEPidRatioEffCont() {
  58   // Destructor
  59
  60   for (Int_t ii = 0; ii < 2; ++ii) {
  61     for (Int_t kk = 0; kk < 4; ++kk)
  62       if (fLabelsRec[ii][kk]) delete[] fLabelsRec[ii][kk];
  63     if (fLabelsRec[ii]) delete[] fLabelsRec[ii];
  64   }
  65   if (fLabelsRec) delete[] fLabelsRec;
  66
  67
  68 }
  69
  70 //________________________________________________________________________
  71 void AliEbyEPidRatioEffCont::Process() {
  72   // -- Process event
  73
  74   // -- Setup (clean, create and fill) MC labels
  75   FillMCLabels();
  76
  77   // -- Fill  MC histograms for efficiency studies
  78   FillMCEffHist();
  79
  80   return;
  81 }
  82
  83 //________________________________________________________________________
  84 void AliEbyEPidRatioEffCont::Init() {
  85   // -- Init eventwise
  86
  87   fLabelsRec = new Int_t**[2];
  88   for (Int_t ii = 0 ; ii < 2; ++ii) {
  89     fLabelsRec[ii] = new Int_t*[4];
  90     for (Int_t kk = 0 ; kk < 4; ++kk)
  91       fLabelsRec[ii][kk] = NULL;
  92   }
  93 }
  94
  95 //________________________________________________________________________
  96 void AliEbyEPidRatioEffCont::CreateHistograms() {
  97   // Copied from NetParticle class
  98   Int_t    binHnEff[10] = { AliEbyEPidRatioHelper::fgkfHistNBinsCent, 4,
  99                             AliEbyEPidRatioHelper::fgkfHistNBinsSign, 2,      2  ,      2  ,
 100                             AliEbyEPidRatioHelper::fgkfHistNBinsEta,
 101                             AliEbyEPidRatioHelper::fgkfHistNBinsRap,
 102                             AliEbyEPidRatioHelper::fgkfHistNBinsPhi,
 103                             AliEbyEPidRatioHelper::fgkfHistNBinsPt};
 104
 105   Double_t minHnEff[10] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0], -0.5,
 106                            AliEbyEPidRatioHelper::fgkfHistRangeSign[0], -0.5, -0.5, -0.5,
 107                            AliEbyEPidRatioHelper::fgkfHistRangeEta[0],
 108                            AliEbyEPidRatioHelper::fgkfHistRangeRap[0],
 109                            AliEbyEPidRatioHelper::fgkfHistRangePhi[0],
 110                            AliEbyEPidRatioHelper::fgkfHistRangePt[0]};
 111
 112
 113   Double_t maxHnEff[10] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[1], 3.5,
 114                            AliEbyEPidRatioHelper::fgkfHistRangeSign[1], 1.5,      1.5,      1.5,
 115                            AliEbyEPidRatioHelper::fgkfHistRangeEta[1],
 116                            AliEbyEPidRatioHelper::fgkfHistRangeRap[1],
 117                            AliEbyEPidRatioHelper::fgkfHistRangePhi[1],
 118                            AliEbyEPidRatioHelper::fgkfHistRangePt[1]};
 119
 120   fHnEffMc    = new THnSparseF("hnEffMc", "cent:pid:SignMC:findable:recStatus:pidStatus:etaMC:yMC:phiMC:ptMC", 10, binHnEff, minHnEff, maxHnEff);
 121   fHnEffRec   = new THnSparseF("hnEffRec", "cent:pid:SignMC:findable:recStatus:pidStatus:etaRec:yRec:phiRec:ptRec", 10, binHnEff, minHnEff, maxHnEff);
 122   //fHnEffRecMc = new THnSparseF("hnEffRecMMc", "cent:pid:SignMC:findable:recStatus:pidStatus:deltaEta:deltaY:deltaPhi:deltaPt" 10, binHnEff, minHnEff, maxHnEff);
 123
 124   fHnEffMc->Sumw2();
 125   fHnEffRec->Sumw2();
 126
 127
 128
 129   fHnEffMc->GetAxis(0)->SetTitle("centrality");                   //  0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins
 130   fHnEffMc->GetAxis(1)->SetTitle("N_{ch}|N_{#pi}|N_{K}|N_{p}");                //  0 | 1 | 2 | 3
 131   fHnEffMc->GetAxis(2)->SetTitle("sign");                         //  -1 | 0 | +1
 132   fHnEffMc->GetAxis(3)->SetTitle("findable");                     //  0 not findable      |  1 findable
 133   fHnEffMc->GetAxis(4)->SetTitle("recStatus");                    //  0 not reconstructed |  1 reconstructed
 134   fHnEffMc->GetAxis(5)->SetTitle("recPid");                       //  0 not accepted      |  1 accepted
 135   fHnEffMc->GetAxis(6)->SetTitle("#eta_{MC}");                    //  eta  [-0.9, 0.9]
 136   fHnEffMc->GetAxis(7)->SetTitle("#it{y}_{MC}");                  //  rapidity  [-0.5, 0.5]
 137   fHnEffMc->GetAxis(8)->SetTitle("#varphi_{MC} (rad)");           //  phi  [ 0. , 2Pi]
 138   fHnEffMc->GetAxis(9)->SetTitle("#it{p}_{T,MC} (GeV/#it{c})");   //  pT   [ 0.2, 2.3]
 139
 140
 141   fHnEffRec->GetAxis(0)->SetTitle("centrality");                   //  0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins
 142   fHnEffRec->GetAxis(1)->SetTitle("N_{ch}|N_{#pi}|N_{K}|N_{p}");                //  0 | 1 | 2 | 3
 143   fHnEffRec->GetAxis(2)->SetTitle("sign");                         //  -1 | 0 | +1
 144   fHnEffRec->GetAxis(3)->SetTitle("findable");                     //  0 not findable      |  1 findable
 145   fHnEffRec->GetAxis(4)->SetTitle("recStatus");                    //  0 not reconstructed |  1 reconstructed
 146   fHnEffRec->GetAxis(5)->SetTitle("recPid");                       //  0 not accepted      |  1 accepted
 147   fHnEffRec->GetAxis(6)->SetTitle("#eta_{Rec}");                   //  eta  [-0.9, 0.9]
 148   fHnEffRec->GetAxis(7)->SetTitle("#it{y}_{Rec}");                 //  rapidity  [-0.5, 0.5]
 149   fHnEffRec->GetAxis(8)->SetTitle("#varphi_{Rec} (rad)");          //  phi  [ 0. , 2Pi]
 150   fHnEffRec->GetAxis(9)->SetTitle("#it{p}_{T,Rec} (GeV/#it{c})");  //  pt   [ 0.2, 2.3]
 151
 152   fHelper->BinLogAxis(fHnEffMc, 9);
 153   fHelper->BinLogAxis(fHnEffRec, 9);
 154
 155   /*
 156      >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Copied from NetParticle
 157
 158      creation -> findable -> reconstructed -> pid_TPC+TOF
 159         (1)         (2)          (3)            (4)
 160                                                                       ||   findable | recStatus | recPid
 161      1) all primary probeParticles_MC                                 ||      -           -         -
 162      2) all findable primary probeParticles_MC                        ||      x           -         -
 163      3) all reconstructed primary probeParticles_MC                   ||      x           x         -
 164      4) all reconstructed primary probeParticles_MC & recPid_TPC+TOF  ||      x           x         x
 165
 166      <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 167   */
 168
 169   // ------------------------------------------------------------------
 170   // -- Create THnSparse - Cont
 171   // ------------------------------------------------------------------
 172
 173   Int_t    binHnCont[8] = {AliEbyEPidRatioHelper::fgkfHistNBinsCent, 4,
 174                            AliEbyEPidRatioHelper::fgkfHistNBinsSign, 8,
 175                            AliEbyEPidRatioHelper::fgkfHistNBinsEta,
 176                            AliEbyEPidRatioHelper::fgkfHistNBinsRap,
 177                            AliEbyEPidRatioHelper::fgkfHistNBinsPhi,
 178                            AliEbyEPidRatioHelper::fgkfHistNBinsPt};
 179
 180   Double_t minHnCont[8] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0], -0.5,
 181                            AliEbyEPidRatioHelper::fgkfHistRangeSign[0],  0.5,
 182                            AliEbyEPidRatioHelper::fgkfHistRangeEta[0],
 183                            AliEbyEPidRatioHelper::fgkfHistRangeRap[0],
 184                            AliEbyEPidRatioHelper::fgkfHistRangePhi[0],
 185                            AliEbyEPidRatioHelper::fgkfHistRangePt[0]};
 186
 187   Double_t maxHnCont[8] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[1], 3.5,
 188                            AliEbyEPidRatioHelper::fgkfHistRangeSign[1], 8.5,
 189                            AliEbyEPidRatioHelper::fgkfHistRangeEta[1],
 190                            AliEbyEPidRatioHelper::fgkfHistRangeRap[1],
 191                            AliEbyEPidRatioHelper::fgkfHistRangePhi[1],
 192                            AliEbyEPidRatioHelper::fgkfHistRangePt[1]};
 193
 194   fHnContMc  = new THnSparseF("hnContMc", "cent:pid:SignMC:contPart:etaMC:yMC:phiMC:ptMC",8, binHnCont, minHnCont, maxHnCont);
 195   fHnContRec = new THnSparseF("hnContRec", "cent:pid:SignRec:contPart:etaRec:yRec:phiRec:ptRec",8, binHnCont, minHnCont, maxHnCont);
 196
 197   fHnContMc->Sumw2();
 198   fHnContRec->Sumw2();
 199
 200   fHnContMc->GetAxis(0)->SetTitle("centrality");                   //  0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins
 201   fHnContMc->GetAxis(1)->SetTitle("N_{ch}|N_{#pi}|N_{K}|N_{p}");                //  0 | 1 | 2 | 3
 202   fHnContMc->GetAxis(2)->SetTitle("sign");                         //  -1 | 0 | +1
 203   fHnContMc->GetAxis(3)->SetTitle("contPart");                     //  1 pi | 2 K | 3 p | 4 e | 5 mu | 6 other | 7 p from WeakDecay | 8 p from Material
 204   fHnContMc->GetAxis(4)->SetTitle("#eta_{MC}");                    //  eta  [-0.9,0.9]
 205   fHnContMc->GetAxis(5)->SetTitle("#it{y}_{MC}");                  //  rapidity  [-0.5, 0.5]
 206   fHnContMc->GetAxis(6)->SetTitle("#varphi_{MC} (rad)");           //  phi  [ 0. ,2Pi]
 207   fHnContMc->GetAxis(7)->SetTitle("#it{p}_{T,MC} (GeV/#it{c})");   //  pT   [ 0.2,2.3]
 208
 209
 210   fHnContRec->GetAxis(0)->SetTitle("centrality");                   //  0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins
 211   fHnContRec->GetAxis(1)->SetTitle("N_{ch}|N_{#pi}|N_{K}|N_{p}");                //  0 | 1 | 2 | 3
 212   fHnContRec->GetAxis(2)->SetTitle("sign");                         //  -1 | 0 | +1
 213   fHnContRec->GetAxis(3)->SetTitle("contPart");                     //  1 pi | 2 K | 3 p | 4 e | 5 mu | 6 other | 7 p from WeakDecay | 8 p from Material
 214   fHnContRec->GetAxis(4)->SetTitle("#eta_{Rec}");                   //  eta  [-0.9, 0.9]
 215   fHnContRec->GetAxis(5)->SetTitle("#it{y}_{Rec}");                 //  rapidity  [-0.5, 0.5]
 216   fHnContRec->GetAxis(6)->SetTitle("#varphi_{Rec} (rad)");          //  phi  [ 0. , 2Pi]
 217   fHnContRec->GetAxis(7)->SetTitle("#it{p}_{T,Rec} (GeV/#it{c})"); //  pt   [ 0.2, 2.3]
 218
 219
 220   //  fHnCont->GetAxis(12)->SetTitle("#eta_{MC}-#eta_{Rec}");                      //  eta  [-0.9, 0.9]
 221   // fHnCont->GetAxis(13)->SetTitle("#it{y}_{MC}-#it{y}_{Rec}");                  //  rapidity  [-0.5, 0.5]
 222   // fHnCont->GetAxis(14)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)");          //  phi  [ -2Pi , 2Pi]
 223   // fHnCont->GetAxis(15)->SetTitle("#it{p}_{T,MC}-#it{p}_{T,Rec} (GeV/#it{c})"); //  pt   [ -2.3, 2.3]
 224   // fHnCont->GetAxis(16)->SetTitle("sign_{MC}-sign_{Rec}");                      //  -2 | 0 | +2
 225   // fHnCont->GetAxis(17)->SetTitle("N_{ch}|N_{#pi}|N_{K}|N_{p}");                //  0 | 1 | 2 | 3
 226
 227   fHelper->BinLogAxis(fHnContMc,  7);
 228   fHelper->BinLogAxis(fHnContRec, 7);
 229
 230   return;
 231 }
 232
 233 //________________________________________________________________________
 234 Int_t AliEbyEPidRatioEffCont::Setup() {
 235   // -- Setup eventwise
 236
 237   // -- Create label arrays
 238   for(Int_t i = 0; i < 4; i++) {
 239   fLabelsRec[0][i] = new Int_t[fNTracks];
 240   if(!fLabelsRec[0][i]) {
 241     AliError("Cannot create fLabelsRec[0]");
 242     return -1;
 243   }
 244
 245   fLabelsRec[1][i] = new Int_t[fNTracks];
 246   if(!fLabelsRec[1][i]) {
 247     AliError("Cannot create fLabelsRec[1] for PID");
 248     return -1;
 249   }
 250
 251   for(Int_t ii = 0; ii < fNTracks; ++ii) {
 252     fLabelsRec[0][i][ii] = 0;
 253     fLabelsRec[1][i][ii] = 0;
 254   }
 255   }
 256   return 0;
 257 }
 258
 259 //________________________________________________________________________
 260 void AliEbyEPidRatioEffCont::Reset() {
 261   // -- Reset eventwise
 262   for(Int_t i = 0; i < 4; i++) {
 263     for (Int_t ii = 0; ii < 2 ; ++ii) {
 264       if (fLabelsRec[ii][i])
 265         delete[] fLabelsRec[ii][i];
 266       fLabelsRec[ii][i] = NULL;
 267     }
 268   }
 269 }
 270
 271 //________________________________________________________________________
 272 void AliEbyEPidRatioEffCont::FillMCLabels() {
 273   Float_t etaRange[2];
 274   fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
 275
 276   Float_t ptRange[2];
 277   fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
 278
 279   // -- Track Loop
 280   for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
 281
 282     AliVTrack *track = (fESD) ? static_cast<AliVTrack*>(fESD->GetTrack(idxTrack)) : static_cast<AliVTrack*>(fAOD->GetTrack(idxTrack));
 283
 284     // -- Check if track is accepted for basic parameters
 285     if (!fHelper->IsTrackAcceptedBasicCharged(track))
 286       continue;
 287
 288     // -- Check if accepted - ESD
 289     if (fESD && !fESDTrackCuts->AcceptTrack(dynamic_cast<AliESDtrack*>(track)))
 290       continue;
 291
 292     // -- Check if accepted - AOD
 293     if (fAOD){
 294       AliAODTrack * trackAOD = dynamic_cast<AliAODTrack*>(track);
 295
 296       if (!trackAOD) {
 297         AliError("Pointer to dynamic_cast<AliAODTrack*>(track) = ZERO");
 298         continue;
 299       }
 300       if (!trackAOD->TestFilterBit(fAODtrackCutBit))
 301         continue;
 302
 303       // -- Check if in pT and eta range (is done in ESDTrackCuts for ESDs)
 304       if(!(track->Pt() > ptRange[0] && track->Pt() <= ptRange[1] && TMath::Abs(track->Eta()) <= etaRange[1]))
 305         continue;
 306     }
 307
 308     Int_t gPdgCode = 0;
 309
 310     Int_t iPid = 0;
 311     Double_t pid[3];
 312     if      (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kPion)))  {  iPid = 1; gPdgCode = 211;}
 313     else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kKaon)))  {  iPid = 2; gPdgCode = 321;}
 314     else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kProton))){  iPid = 3; gPdgCode = 2212;}
 315     else iPid = 0;
 316
 317     //  cout << " --- EFF ---- " << iPid << "  " << gPdgCode << endl;
 318
 319     Double_t yP;
 320     if (!fHelper->IsTrackAcceptedRapidity(track, yP, iPid))
 321       continue;
 322
 323     if (!fHelper->IsTrackAcceptedDCA(track))
 324       continue;
 325
 326     Int_t label  = TMath::Abs(track->GetLabel());
 327
 328     // -- Fill Label of all reconstructed
 329     if(iPid != 0) fLabelsRec[0][0][idxTrack]    = label;
 330     fLabelsRec[0][iPid][idxTrack] = label;
 331
 332     // -- Fill Label of all reconstructed && recPid_TPC+TOF
 333     if(iPid != 0) fLabelsRec[1][0][idxTrack]    = label;
 334     fLabelsRec[1][iPid][idxTrack] = label;
 335
 336     // -- Check for contamination and fill contamination THnSparse
 337     CheckContTrack(track, iPid, gPdgCode);
 338     if(iPid != 0) CheckContTrack(track, 0, 0);
 339
 340   } // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
 341
 342   return;
 343 }
 344
 345 //________________________________________________________________________
 346 void AliEbyEPidRatioEffCont::CheckContTrack(AliVTrack *track, Int_t iPid, Int_t gPdgCode) {
 347   Int_t label     = TMath::Abs(track->GetLabel());
 348   Float_t signRec = track->Charge();
 349
 350
 351   AliVParticle* particle = (fESD) ? fMCEvent->GetTrack(label) : static_cast<AliVParticle*>(fArrayMC->At(label));
 352   if (!particle)
 353     return;
 354
 355   Bool_t isPhysicalPrimary = (fESD) ? fStack->IsPhysicalPrimary(label): (static_cast<AliAODMCParticle*>(particle))->IsPhysicalPrimary();
 356   if (iPid == 0) {
 357     if (particle->PdgCode() == (signRec*gPdgCode))
 358       if (isPhysicalPrimary)
 359         return;
 360   }
 361   else {
 362     if (isPhysicalPrimary)
 363       return;
 364   }
 365
 366   // -- Check if secondaries from material or weak decay
 367   Bool_t isSecondaryFromWeakDecay = (fESD) ? fStack->IsSecondaryFromWeakDecay(label) : (static_cast<AliAODMCParticle*>(particle))->IsSecondaryFromWeakDecay();
 368   Bool_t isSecondaryFromMaterial  = (fESD) ? fStack->IsSecondaryFromMaterial(label)  : (static_cast<AliAODMCParticle*>(particle))->IsSecondaryFromMaterial();
 369
 370   // -- Get PDG Charge of contaminating particle
 371   Float_t signMC = 0.;
 372   if      (particle->Charge() == 0.) signMC =  0.;
 373   else if (particle->Charge() <  0.) signMC = -1.;
 374   else if (particle->Charge() >  0.) signMC =  1.;
 375
 376   // -- Get contaminating particle
 377   Float_t contPart = 0;
 378   if        (isSecondaryFromWeakDecay)                contPart = 7; // probeParticle from WeakDecay
 379   else if   (isSecondaryFromMaterial)                 contPart = 8; // probeParticle from Material
 380   else {
 381     if      (TMath::Abs(particle->PdgCode()) ==  211) contPart = 1; // pion
 382     else if (TMath::Abs(particle->PdgCode()) ==  321) contPart = 2; // kaon
 383     else if (TMath::Abs(particle->PdgCode()) == 2212) contPart = 3; // proton
 384     else if (TMath::Abs(particle->PdgCode()) ==   11) contPart = 4; // electron
 385     else if (TMath::Abs(particle->PdgCode()) ==   13) contPart = 5; // muon
 386     else                                              contPart = 6; // other
 387   }
 388
 389   // cout << " --- CONT ---- " << iPid << "  " << gPdgCode << endl;
 390
 391   // -- Get Reconstructed y
 392   //    yRec = y for identified particles | yRec = eta for charged particles
 393   Double_t yRec  = 0.;
 394   fHelper->IsTrackAcceptedRapidity(track, yRec, iPid);
 395
 396   Double_t deltaPhi = particle->Phi()-track->Phi();
 397   if (TMath::Abs(deltaPhi) > TMath::TwoPi()) {
 398     if (deltaPhi < 0)
 399       deltaPhi += TMath::TwoPi();
 400     else
 401       deltaPhi -= TMath::TwoPi();
 402   }
 403
 404   Double_t hnContMc[8]  = {fCentralityBin,iPid,signMC,contPart,particle->Eta(),particle->Y(),particle->Phi(),particle->Pt()};
 405   Double_t hnContRec[8] = {fCentralityBin,iPid,signRec,contPart, track->Eta(),yRec,track->Phi(),track->Pt()};
 406   fHnContMc->Fill(hnContMc);
 407   fHnContRec->Fill(hnContRec);
 408
 409 }
 410
 411 //________________________________________________________________________
 412 void AliEbyEPidRatioEffCont::FillMCEffHist() {
 413   // Fill efficiency THnSparse for ESDs
 414
 415   Float_t etaRange[2];
 416   fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
 417
 418   Int_t nPart  = (fESD) ? fStack->GetNprimary() : fArrayMC->GetEntriesFast();
 419
 420   for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
 421     AliVParticle* particle = (fESD) ? fMCEvent->GetTrack(idxMC) : static_cast<AliVParticle*>(fArrayMC->At(idxMC));
 422
 423     // -- Check basic MC properties -> charged physical primary
 424     if (!fHelper->IsParticleAcceptedBasicCharged(particle, idxMC))
 425       continue;
 426
 427     Int_t iPid = 0;
 428     Int_t gPdgCode = 0;
 429     if ( TMath::Abs(particle->PdgCode())      ==  211 ) {  iPid = 1; gPdgCode = 211;}
 430     else if ( TMath::Abs(particle->PdgCode()) ==  321 ) {  iPid = 2; gPdgCode = 321;}
 431     else if ( TMath::Abs(particle->PdgCode()) == 2212 ) {  iPid = 3; gPdgCode = 2212;}
 432     else {iPid = 0; gPdgCode = 0;}
 433
 434     // -- Check if accepted in rapidity window -- for identified particles
 435     Double_t yMC;
 436     if (iPid != 0) {
 437       if(!fHelper->IsParticleAcceptedRapidity(particle, yMC, iPid))
 438         continue;
 439     } else {
 440       // -- Check if accepted in eta window -- for charged particles
 441       if (TMath::Abs(particle->Eta()) > etaRange[1])
 442         continue;
 443     }
 444
 445     // cout << particle->PdgCode() << "  " <<iPid << endl;
 446
 447     // -- Check if probeParticle / anti-probeParticle
 448     //    > skip check if PID is not required
 449     if (iPid != 0) {
 450       if (TMath::Abs(particle->PdgCode()) != gPdgCode)
 451       continue;
 452     }
 453
 454     // -- Get sign of particle
 455     Float_t signMC    = (particle->PdgCode() < 0) ? -1. : 1.;
 456
 457     // -- Get if particle is findable --- not availible for AODs yet
 458     Float_t findable  = (fESD) ? Float_t(fHelper->IsParticleFindable(idxMC)) : 1.;
 459
 460     // cout << findable << "  " << fHelper->IsParticleFindable(idxMC)<< endl;
 461
 462     // -- Get recStatus and pidStatus
 463     Float_t recStatus = 0.;
 464     Float_t recPid    = 0.;
 465
 466     // -- Get Reconstructed values
 467     Float_t etaRec  = 0.;
 468     Float_t phiRec  = 0.;
 469     Float_t ptRec   = 0.;
 470     Double_t yRec   = 0.;
 471     Float_t signRec = 0.;
 472
 473     // -- Loop over all labels
 474     for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {
 475       if (idxMC == fLabelsRec[0][iPid][idxRec]) {
 476         recStatus = 1.;
 477
 478         if (idxMC == fLabelsRec[1][iPid][idxRec])
 479           recPid = 1.;
 480
 481         AliVTrack *track = NULL;
 482         if(fESD)
 483           track = fESD->GetTrack(idxRec);
 484         else if(fAOD)
 485           track = fAOD->GetTrack(idxRec);
 486
 487         if (track) {
 488           // if no track present (which should not happen)
 489           // -> pt = 0. , which is not in the looked at range
 490
 491           // -- Get Reconstructed values
 492           etaRec  = track->Eta();
 493           phiRec  = track->Phi();
 494           ptRec   = track->Pt();
 495           signRec = track->Charge();
 496           fHelper->IsTrackAcceptedRapidity(track, yRec, iPid); // yRec = y for identified particles | yRec = eta for charged particles
 497         }
 498         break;
 499       }
 500     } // for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {
 501
 502     Double_t deltaPhi = particle->Phi()-phiRec;
 503     if (TMath::Abs(deltaPhi) > TMath::TwoPi()) {
 504       if (deltaPhi < 0)
 505         deltaPhi += TMath::TwoPi();
 506       else
 507         deltaPhi -= TMath::TwoPi();
 508     }
 509
 510     // if (signRec == 0) continue;
 511
 512     if(iPid != 0) {
 513       Double_t hnEffMc[10]  = {fCentralityBin,0,signMC,findable, recStatus,recPid,particle->Eta(), particle->Y(), particle->Phi(),particle->Pt()};
 514       Double_t hnEffRec[10] = {fCentralityBin,0,signRec,findable, recStatus,recPid,etaRec, yRec, phiRec, ptRec};
 515       fHnEffMc->Fill(hnEffMc);
 516       fHnEffRec->Fill(hnEffRec);
 517     }
 518     Double_t hnEffMc[10]  = {fCentralityBin,iPid,signMC,findable, recStatus,recPid,particle->Eta(), particle->Y(), particle->Phi(),particle->Pt()};
 519     Double_t hnEffRec[10] = {fCentralityBin,iPid,signRec,findable, recStatus,recPid,etaRec, yRec, phiRec, ptRec};
 520     fHnEffMc->Fill(hnEffMc);
 521     fHnEffRec->Fill(hnEffRec);
 522
 523
 524     //  cout << signMC << "  " << signRec << "  " << iPid << "  " << gPdgCode << endl;
 525
 526
 527   } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
 528
 529   return;
 530 }