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