PWGCF/EBYE/NetParticle/AliAnalysisNetParticleEffCont.cxx

   1 //-*- Mode: C++ -*-
   2
   3 #include "TMath.h"
   4 #include "TAxis.h"
   5 #include "TDatabasePDG.h"
   6
   7 #include "AliESDEvent.h"
   8 #include "AliESDInputHandler.h"
   9 #include "AliStack.h"
  10 #include "AliMCEvent.h"
  11 #include "AliESDtrackCuts.h"
  12 #include "AliAODEvent.h"
  13 #include "AliAODInputHandler.h"
  14 #include "AliAODMCParticle.h"
  15
  16 #include "AliAnalysisNetParticleEffCont.h"
  17
  18 using namespace std;
  19
  20 // Task for NetParticle checks
  21 // Author: Jochen Thaeder <jochen@thaeder.de>
  22
  23 ClassImp(AliAnalysisNetParticleEffCont)
  24
  25 /*
  26  * ---------------------------------------------------------------------------------
  27  *                            Constructor / Destructor
  28  * ---------------------------------------------------------------------------------
  29  */
  30
  31 //________________________________________________________________________
  32 AliAnalysisNetParticleEffCont::AliAnalysisNetParticleEffCont() :
  33   fHelper(NULL),
  34   fPdgCode(2212),
  35
  36   fESD(NULL),
  37   fESDTrackCuts(NULL),
  38
  39   fAOD(NULL),
  40   fArrayMC(NULL),
  41
  42   fCentralityBin(-1.),
  43   fNTracks(0),
  44
  45   fAODtrackCutBit(1024),
  46
  47   fStack(NULL),
  48   fMCEvent(NULL),
  49
  50   fLabelsRec(NULL),
  51
  52   fHnEff(NULL),
  53   fHnCont(NULL) {
  54   // Constructor
  55
  56   AliLog::SetClassDebugLevel("AliAnalysisNetParticleEffCont",10);
  57 }
  58
  59 //________________________________________________________________________
  60 AliAnalysisNetParticleEffCont::~AliAnalysisNetParticleEffCont() {
  61   // Destructor
  62
  63   if (fLabelsRec){
  64
  65     if (fLabelsRec[0])
  66       delete[] (fLabelsRec[0]);
  67
  68     if (fLabelsRec[1])
  69       delete[] (fLabelsRec[1]);
  70
  71     delete[] fLabelsRec;
  72   }
  73 }
  74
  75 /*
  76  * ---------------------------------------------------------------------------------
  77  *                                 Public Methods
  78  * ---------------------------------------------------------------------------------
  79  */
  80
  81 //________________________________________________________________________
  82 void AliAnalysisNetParticleEffCont::Initialize(AliESDtrackCuts *cuts , AliAnalysisNetParticleHelper* helper, Int_t trackCutBit) {
  83
  84   // -- ESD track cuts
  85   // -------------------
  86   fESDTrackCuts     = cuts;
  87
  88   // -- Get Helper
  89   // ---------------
  90   fHelper           = helper;
  91
  92   // -- Get particle species / pdgCode
  93   // -------------------------
  94   fPdgCode          = AliPID::ParticleCode(fHelper->GetParticleSpecies());
  95
  96   // -- MC Labels for efficiency
  97   // -----------------------------
  98   fLabelsRec        = new Int_t*[2];
  99   for (Int_t ii = 0; ii < 2 ; ++ii)
 100     fLabelsRec[ii] = NULL;
 101
 102   // -- AOD track filter bit
 103   // -------------------------
 104   fAODtrackCutBit = trackCutBit;
 105
 106   // -- Create THnSparse Histograms
 107   // --------------------------------
 108   CreateHistograms();
 109
 110   return;
 111 }
 112
 113 /*
 114  * ---------------------------------------------------------------------------------
 115  *                            Setup/Reset Methods - private
 116  * ---------------------------------------------------------------------------------
 117  */
 118
 119 //________________________________________________________________________
 120 Int_t AliAnalysisNetParticleEffCont::SetupEvent(AliESDInputHandler *esdHandler, AliAODInputHandler *aodHandler, AliMCEvent *mcEvent) {
 121   // -- Setup event
 122
 123   // -- Reset of event
 124   // -------------------
 125   ResetEvent();
 126
 127   // -- Setup of event
 128   // -------------------
 129
 130   // -- Get ESD objects
 131   if(esdHandler) {
 132     fESD     = esdHandler->GetEvent();
 133     fNTracks = fESD->GetNumberOfTracks();
 134   }
 135
 136   // -- Get AOD objects
 137   else if(aodHandler) {
 138     fAOD     = aodHandler->GetEvent();
 139     fNTracks = fAOD->GetNumberOfTracks();
 140
 141     fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
 142     if (!fArrayMC)
 143       AliFatal("No array of MC particles found !!!"); // MW  no AliFatal use return values
 144   }
 145
 146   // -- Get MC objects
 147   if (mcEvent) {
 148     fMCEvent     = mcEvent;
 149     if (fMCEvent)
 150       fStack     = fMCEvent->Stack();
 151   }
 152
 153   fCentralityBin = fHelper->GetCentralityBin();
 154
 155   // -- Create label arrays
 156   // ------------------------
 157   fLabelsRec[0] = new Int_t[fNTracks];
 158   if(!fLabelsRec[0]) {
 159     AliError("Cannot create fLabelsRec[0]");
 160     return -1;
 161   }
 162
 163   fLabelsRec[1] = new Int_t[fNTracks];
 164   if(!fLabelsRec[1]) {
 165     AliError("Cannot create fLabelsRec[1] for TPC");
 166     return -1;
 167  }
 168
 169   for(Int_t ii = 0; ii < fNTracks; ++ii) {
 170     fLabelsRec[0][ii] = 0;
 171     fLabelsRec[1][ii] = 0;
 172   }
 173
 174   return 0;
 175 }
 176
 177 //________________________________________________________________________
 178 void AliAnalysisNetParticleEffCont::ResetEvent() {
 179   // -- Reset event
 180
 181   for (Int_t ii = 0; ii < 2 ; ++ii) {
 182     if (fLabelsRec[ii])
 183       delete[] fLabelsRec[ii];
 184     fLabelsRec[ii] = NULL;
 185   }
 186
 187   return;
 188 }
 189
 190 //________________________________________________________________________
 191 void AliAnalysisNetParticleEffCont::Process() {
 192   // -- Process event
 193
 194   // -- Setup (clean, create and fill) MC labels
 195   // ---------------------------------------------
 196   FillMCLabels();
 197
 198   // -- Fill  MC histograms for efficiency studies
 199   // -----------------------------------------------
 200   if(fESD)
 201     FillMCEffHist();
 202   else if(fAOD)
 203     FillMCEffHistAOD();
 204
 205   return;
 206 }
 207
 208 /*
 209  * ---------------------------------------------------------------------------------
 210  *                                 Private Methods
 211  * ---------------------------------------------------------------------------------
 212  */
 213
 214 //________________________________________________________________________
 215 void AliAnalysisNetParticleEffCont::CreateHistograms() {
 216   // -- Create histograms
 217
 218   Double_t dCent[2] = {-0.5, 8.5};
 219   Int_t iCent       = 9;
 220
 221   Double_t dEta[2]  = {-0.9, 0.9};
 222   Int_t iEta        = Int_t((dEta[1]-dEta[0]) / 0.075) +1 ;
 223
 224   Double_t dRap[2]  = {-0.5, 0.5};
 225   Int_t iRap        = Int_t((dRap[1]-dRap[0]) / 0.075) +1 ;
 226
 227   Double_t dPhi[2]  = {0.0, TMath::TwoPi()};
 228   Int_t iPhi        = 42;
 229
 230   Double_t dPt[2]   = {0.1, 3.0};
 231   Int_t iPt         = Int_t((dPt[1]-dPt[0]) / 0.075);
 232
 233   Double_t dSign[2] = {-1.5, 1.5};
 234   Int_t iSign       = 3;
 235
 236   // ------------------------------------------------------------------
 237   // -- Create THnSparseF - Eff
 238   // ------------------------------------------------------------------
 239
 240   //                           cent:   etaMC:     yMC:  phiMC:   ptMC:     sign:findable:recStatus:pidStatus:   etaRec:  phiRec:  ptRec:deltaEta: deltaPhi: deltaPt
 241   Int_t    binHnEff[15] = {   iCent,    iEta,    iRap,    iPhi,    iPt,    iSign,       2,      2  ,      2  ,    iEta,    iPhi,    iPt,    iEta, 2*iPhi+1, 2*iPt+1};
 242   Double_t minHnEff[15] = {dCent[0], dEta[0], dRap[0], dPhi[0], dPt[0], dSign[0],    -0.5,     -0.5,     -0.5, dEta[0], dPhi[0], dPt[0], dEta[0], -dPhi[1], -dPt[1]};
 243   Double_t maxHnEff[15] = {dCent[1], dEta[1], dRap[1], dPhi[1], dPt[1], dSign[1],     1.5,      1.5,      1.5, dEta[1], dPhi[1], dPt[1], dEta[1],  dPhi[1],  dPt[1]};
 244
 245   fHnEff = new THnSparseF("fHnEff", "cent:etaMC:yMC:phiMC:ptMC:sign:findable:recStatus:pidStatus:etaRec:phiRec:ptRec:deltaEta:deltaPhi:deltaPt",
 246                           15, binHnEff, minHnEff, maxHnEff);
 247
 248   fHnEff->Sumw2();
 249
 250   fHnEff->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
 251   fHnEff->GetAxis(1)->SetTitle("#eta_{MC}");                    //  eta  [-0.9, 0.9]
 252   fHnEff->GetAxis(2)->SetTitle("#it{y}_{MC}");                  //  rapidity  [-0.5, 0.5]
 253   fHnEff->GetAxis(3)->SetTitle("#varphi_{MC} (rad)");           //  phi  [ 0. , 2Pi]
 254   fHnEff->GetAxis(4)->SetTitle("#it{p}_{T,MC} (GeV/#it{c})");   //  pt   [ 0.1, 3.0]
 255
 256   fHnEff->GetAxis(5)->SetTitle("sign");                         //  -1 | 0 | +1
 257   fHnEff->GetAxis(6)->SetTitle("findable");                     //  0 not findable      |  1 findable
 258   fHnEff->GetAxis(7)->SetTitle("recStatus");                    //  0 not reconstructed |  1 reconstructed
 259   fHnEff->GetAxis(8)->SetTitle("recPid");                       //  0 not accepted      |  1 accepted
 260
 261   fHnEff->GetAxis(9)->SetTitle("#eta_{Rec}");                   //  eta  [-0.9, 0.9]
 262   fHnEff->GetAxis(10)->SetTitle("#varphi_{Rec} (rad)");         //  phi  [ 0. , 2Pi]
 263   fHnEff->GetAxis(11)->SetTitle("#it{p}_{T,Rec} (GeV/#it{c})"); //  pt   [ 0.1, 3.0]
 264
 265   fHnEff->GetAxis(12)->SetTitle("#eta_{MC}-#eta_{Rec}");                      //  eta  [-0.9,0.9]
 266   fHnEff->GetAxis(13)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)");          //  phi  [ 0. ,2Pi]
 267   fHnEff->GetAxis(14)->SetTitle("#it{p}_{T,MC}-#it{p}_{T,Rec} (GeV/#it{c})"); //  pt   [ 3.0,3.0]
 268
 269   fHelper->BinLogAxis(fHnEff,  4);
 270   fHelper->BinLogAxis(fHnEff, 11);
 271
 272   /*
 273      >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 274
 275      creation -> findable -> reconstructed -> pid_TPC+TOF
 276         (1)         (2)          (3)            (4)
 277                                                                       ||   findable | recStatus | recPid
 278      1) all primary probeParticles_MC                                 ||      -           -         -
 279      2) all findable primary probeParticles_MC                        ||      x           -         -
 280      3) all reconstructed primary probeParticles_MC                   ||      x           x         -
 281      4) all reconstructed primary probeParticles_MC & recPid_TPC+TOF  ||      x           x         x
 282
 283      <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 284   */
 285
 286   // ------------------------------------------------------------------
 287   // -- Create THnSparseF - Cont
 288   // ------------------------------------------------------------------
 289
 290   //                            cent:   etaMC:    yMC:   phiMC:   ptMC:     sign:contPart: contSign:  etaRec:  phiRec:  ptRec:deltaEta: deltaPhi: deltaPt
 291   Int_t    binHnCont[14] = {   iCent,    iEta,    iRap,    iPhi,    iPt,    iSign,       8,    iSign,    iEta,    iPhi,    iPt,    iEta, 2*iPhi+1, 2*iPt+1};
 292   Double_t minHnCont[14] = {dCent[0], dEta[0], dRap[0], dPhi[0], dPt[0], dSign[0],     0.5, dSign[0], dEta[0], dPhi[0], dPt[0], dEta[0], -dPhi[1], -dPt[1]};
 293   Double_t maxHnCont[14] = {dCent[1], dEta[1], dRap[1], dPhi[1], dPt[1], dSign[1],     8.5, dSign[1], dEta[1], dPhi[1], dPt[1], dEta[1],  dPhi[1],  dPt[1]};
 294   fHnCont = new THnSparseF("fHnCont", "cent:etaMC:yMC:phiMC:ptMC:sign:contPart:contSign:etaRec:phiRec:ptRec:deltaEta:deltaPhi:deltaPt",
 295                            14, binHnCont, minHnCont, maxHnCont);
 296
 297   fHnCont->Sumw2();
 298
 299   fHnCont->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
 300   fHnCont->GetAxis(1)->SetTitle("#eta_{MC}");                    //  eta  [-0.9,0.9]
 301   fHnCont->GetAxis(2)->SetTitle("#it{y}_{MC}");                  //  rapidity  [-0.5, 0.5]
 302   fHnCont->GetAxis(3)->SetTitle("#varphi_{MC} (rad)");           //  phi  [ 0. ,2Pi]
 303   fHnCont->GetAxis(4)->SetTitle("#it{p}_{T,MC} (GeV/#it{c})");   //  pT   [ 0.1,1.3]
 304   fHnCont->GetAxis(5)->SetTitle("sign");                         //  -1 | 0 | +1
 305
 306   fHnCont->GetAxis(6)->SetTitle("contPart");                     //  1 pi | 2 K | 3 p | 4 e | 5 mu | 6 other | 7 p from WeakDecay | 8 p from Material
 307   fHnCont->GetAxis(7)->SetTitle("contSign");                     //  -1 | 0 | +1
 308
 309   fHnCont->GetAxis(8)->SetTitle("#eta_{Rec}");                   //  eta  [-0.9, 0.9]
 310   fHnCont->GetAxis(9)->SetTitle("#varphi_{Rec} (rad)");          //  phi  [ 0. , 2Pi]
 311   fHnCont->GetAxis(10)->SetTitle("#it{p}_{T,Rec} (GeV/#it{c})"); //  pt   [ 0.1, 3.0]
 312
 313   fHnCont->GetAxis(11)->SetTitle("#eta_{MC}-#eta_{Rec}");                      //  eta  [-0.9,0.9]
 314   fHnCont->GetAxis(12)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)");          //  phi  [ 0. ,2Pi]
 315   fHnCont->GetAxis(13)->SetTitle("#it{p}_{T,MC}-#it{p}_{T,Rec} (GeV/#it{c})"); //  pt   [ 3.0,3.0]
 316
 317   fHelper->BinLogAxis(fHnCont,  4);
 318   fHelper->BinLogAxis(fHnCont, 10);
 319
 320   // ------------------------------------------------------------------
 321
 322   return;
 323 }
 324
 325 //________________________________________________________________________
 326 void AliAnalysisNetParticleEffCont::FillMCLabels() {
 327   // Fill MC labels
 328   // Loop over ESD tracks and fill arrays with MC lables
 329   //  fLabelsRec[0] : all Tracks
 330   //  fLabelsRec[1] : all Tracks accepted by PID of TPC
 331   // Check every accepted track if correctly identified
 332   //  otherwise check for contamination
 333
 334   for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
 335
 336     AliVTrack *track = (fESD) ? static_cast<AliVTrack*>(fESD->GetTrack(idxTrack)) : static_cast<AliVTrack*>(fAOD->GetTrack(idxTrack));
 337
 338     // -- Check if track is accepted for basic parameters
 339     if (!fHelper->IsTrackAcceptedBasicCharged(track))
 340       continue;
 341
 342     // -- Check if accepted - ESD
 343     if (fESD && !fESDTrackCuts->AcceptTrack(dynamic_cast<AliESDtrack*>(track)))
 344       continue;
 345
 346     // -- Check if accepted - AOD
 347     if (fAOD && (dynamic_cast<AliAODTrack*>(track))){
 348       if(!((dynamic_cast<AliAODTrack*>(track))->TestFilterBit(fAODtrackCutBit)))
 349         continue;
 350     }
 351
 352     // -- Check if accepted in rapidity window
 353     Double_t yP;
 354     if (!fHelper->IsTrackAcceptedRapidity(track, yP))
 355       continue;
 356
 357     // -- Check if accepted with thighter DCA cuts
 358     // -- returns kTRUE for AODs for now : MW
 359     if (!fHelper->IsTrackAcceptedDCA(track))
 360       continue;
 361
 362     Int_t label  = TMath::Abs(track->GetLabel());
 363
 364     // -- Fill Label of all reconstructed
 365     fLabelsRec[0][idxTrack] = label;
 366
 367     // -- Check if accepted by PID from TPC or TPC+TOF
 368     Double_t pid[2];
 369     if (!fHelper->IsTrackAcceptedPID(track, pid))
 370       continue;
 371
 372     // -- Fill Label of all reconstructed && recPid_TPC+TOF
 373     fLabelsRec[1][idxTrack] = label;
 374
 375     // -- Check for contamination and fill contamination THnSparse
 376     if (fESD)
 377       CheckContTrack(label, track->Charge(), idxTrack);
 378     else if (fAOD)
 379       CheckContTrackAOD(label, track->Charge(), idxTrack);
 380
 381   } // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
 382
 383   return;
 384 }
 385
 386 //________________________________________________________________________
 387 void AliAnalysisNetParticleEffCont::CheckContTrack(Int_t label, Float_t sign, Int_t idxTrack) {
 388   // Check if particle is contamination or correctly identified for ESDs
 389   // Fill contamination THnSparse
 390
 391   TParticle* particle = fStack->Particle(label);
 392   if (!particle)
 393     return;
 394
 395   Bool_t isSecondaryFromWeakDecay = kFALSE;
 396   Bool_t isSecondaryFromMaterial  = kFALSE;
 397
 398   // -- Check if correctly identified
 399   //    > Check for Primaries and Secondaries
 400   if (fHelper->GetUsePID()) {
 401     if (particle->GetPdgCode() == (sign*fPdgCode)) {
 402
 403       // Check if is physical primary -> all ok
 404       if (fStack->IsPhysicalPrimary(label))
 405         return;
 406
 407       // -- Check if secondaries from material or weak decay
 408       isSecondaryFromWeakDecay = fStack->IsSecondaryFromWeakDecay(label);
 409       isSecondaryFromMaterial  = fStack->IsSecondaryFromMaterial(label);
 410     }
 411   }
 412
 413   // -- If no PID is required
 414   //    > only check for Primaries and Secondaries
 415   else {
 416     // Check if is physical primary -> all ok
 417     if (fStack->IsPhysicalPrimary(label))
 418       return;
 419
 420     // -- Check if secondaries from material or weak decay
 421     isSecondaryFromWeakDecay = fStack->IsSecondaryFromWeakDecay(label);
 422     isSecondaryFromMaterial  = fStack->IsSecondaryFromMaterial(label);
 423   }
 424
 425   // -- Get PDG Charge
 426   Float_t contSign = 0.;
 427   if      (particle->GetPDG()->Charge() == 0.) contSign =  0.;
 428   else if (particle->GetPDG()->Charge() <  0.) contSign = -1.;
 429   else if (particle->GetPDG()->Charge() >  0.) contSign =  1.;
 430
 431   // -- Get contaminating particle
 432   Float_t contPart = 0;
 433   if        (isSecondaryFromWeakDecay)                   contPart = 7; // probeParticle from WeakDecay
 434   else if   (isSecondaryFromMaterial)                    contPart = 8; // probeParticle from Material
 435   else {
 436     if      (TMath::Abs(particle->GetPdgCode()) ==  211) contPart = 1; // pion
 437     else if (TMath::Abs(particle->GetPdgCode()) ==  321) contPart = 2; // kaon
 438     else if (TMath::Abs(particle->GetPdgCode()) == 2212) contPart = 3; // proton
 439     else if (TMath::Abs(particle->GetPdgCode()) ==   11) contPart = 4; // electron
 440     else if (TMath::Abs(particle->GetPdgCode()) ==   13) contPart = 5; // muon
 441     else                                                 contPart = 6; // other
 442   }
 443
 444   // -- Get Reconstructed values
 445   Float_t etaRec = 0.;
 446   Float_t phiRec = 0.;
 447   Float_t ptRec  = 0.;
 448
 449   AliESDtrack *track = fESD->GetTrack(idxTrack);
 450   if (track) {
 451     // if no track present (which should not happen)
 452     // -> pt = 0. , which is not in the looked at range
 453
 454     // -- Get Reconstructed eta/phi/pt
 455     etaRec = track->Eta();
 456     phiRec = track->Phi();
 457     ptRec  = track->Pt();
 458   }
 459
 460   // -- Fill THnSparse
 461   Double_t hnCont[14] = {fCentralityBin, particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), sign,
 462                          contPart, contSign, etaRec, phiRec, ptRec,
 463                          particle->Eta()-etaRec, particle->Phi()-phiRec, particle->Pt()-ptRec};
 464   fHnCont->Fill(hnCont);
 465 }
 466
 467 //________________________________________________________________________
 468 void AliAnalysisNetParticleEffCont::CheckContTrackAOD(Int_t label, Float_t sign, Int_t idxTrack) {
 469   // Check if particle is contamination or correctly identified for AODs
 470   // Fill contamination THnSparse
 471
 472   AliAODMCParticle* particle = static_cast<AliAODMCParticle*>(fArrayMC->At(label));
 473
 474   if (!particle)
 475     return;
 476
 477   Bool_t isSecondaryFromWeakDecay = kFALSE;
 478   Bool_t isSecondaryFromMaterial  = kFALSE;
 479
 480   // -- Check if correctly identified
 481   //    > Check for Primaries and Secondaries
 482   if (fHelper->GetUsePID()) {
 483     if (particle->GetPdgCode() == (sign*fPdgCode)) {
 484
 485       // Check if is physical primary -> all ok
 486       if (particle->IsPhysicalPrimary())
 487         return;
 488
 489       // -- Check if secondaries from material or weak decay
 490       isSecondaryFromWeakDecay = particle->IsSecondaryFromWeakDecay();
 491       isSecondaryFromMaterial  = particle->IsSecondaryFromMaterial();
 492     }
 493   }
 494
 495   // -- If no PID is required
 496   //    > only check for Primaries and Secondaries
 497   else {
 498     // Check if is physical primary -> all ok
 499     if (particle->IsPhysicalPrimary())
 500       return;
 501
 502     // -- Check if secondaries from material or weak decay
 503     isSecondaryFromWeakDecay = particle->IsSecondaryFromWeakDecay();
 504     isSecondaryFromMaterial  = particle->IsSecondaryFromMaterial();
 505   }
 506
 507   // -- Get PDG Charge
 508   Float_t contSign = 0.;
 509   // Crashing (MW)? --> why do we need the PDG charge?
 510   // if      ((TDatabasePDG::Instance()->GetParticle(particle->PdgCode()))->Charge() == 0.) contSign =  0.;
 511   // else if ((TDatabasePDG::Instance()->GetParticle(particle->PdgCode()))->Charge() <  0.) contSign = -1.;
 512   // else if ((TDatabasePDG::Instance()->GetParticle(particle->PdgCode()))->Charge() >  0.) contSign =  1.;
 513    if      (particle->Charge() == 0.) contSign =  0.;
 514    else if (particle->Charge() <  0.) contSign = -1.;
 515    else if (particle->Charge() >  0.) contSign =  1.;
 516
 517   // -- Get contaminating particle
 518   Float_t contPart = 0;
 519   if        (isSecondaryFromWeakDecay)                   contPart = 7; // probeParticle from WeakDecay
 520   else if   (isSecondaryFromMaterial)                    contPart = 8; // probeParticle from Material
 521   else {
 522     if      (TMath::Abs(particle->GetPdgCode()) ==  211) contPart = 1; // pion
 523     else if (TMath::Abs(particle->GetPdgCode()) ==  321) contPart = 2; // kaon
 524     else if (TMath::Abs(particle->GetPdgCode()) == 2212) contPart = 3; // proton
 525     else if (TMath::Abs(particle->GetPdgCode()) ==   11) contPart = 4; // electron
 526     else if (TMath::Abs(particle->GetPdgCode()) ==   13) contPart = 5; // muon
 527     else                                                 contPart = 6; // other
 528   }
 529
 530   // -- Get Reconstructed values
 531   Float_t etaRec = 0.;
 532   Float_t phiRec = 0.;
 533   Float_t ptRec  = 0.;
 534
 535   AliAODTrack *track = fAOD->GetTrack(idxTrack);
 536
 537   if (track) {
 538     // if no track present (which should not happen)
 539     // -> pt = 0. , which is not in the looked at range
 540
 541     // -- Get Reconstructed eta/phi/pt
 542     etaRec = track->Eta();
 543     phiRec = track->Phi();
 544     ptRec  = track->Pt();
 545   }
 546
 547   // -- Fill THnSparse
 548   Double_t hnCont[14] = {fCentralityBin, particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), sign,
 549                          contPart, contSign, etaRec, phiRec, ptRec,
 550                          particle->Eta()-etaRec, particle->Phi()-phiRec, particle->Pt()-ptRec};
 551   fHnCont->Fill(hnCont);
 552 }
 553
 554 //________________________________________________________________________
 555 void AliAnalysisNetParticleEffCont::FillMCEffHist() {
 556   // Fill efficiency THnSparse for ESDs
 557
 558   Int_t nPart  = fStack->GetNprimary();
 559
 560   Float_t etaRange[2];
 561   fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
 562
 563   for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
 564     TParticle* particle = fStack->Particle(idxMC);
 565
 566     // -- Check basic MC properties -> charged physical primary
 567     if (!fHelper->IsParticleAcceptedBasicCharged(particle, idxMC))
 568       continue;
 569
 570     // -- Check rapidity window
 571     Double_t yMC;
 572     if (!fHelper->GetUsePID())
 573       yMC = etaRange[1];
 574     if (!fHelper->IsParticleAcceptedRapidity(particle, yMC))
 575       continue;
 576
 577     // -- Check if probeParticle / anti-probeParticle
 578     //    > skip check if PID is not required
 579     if (fHelper->GetUsePID() && TMath::Abs(particle->GetPdgCode()) != fPdgCode){
 580       printf("SHOULD NOT HAPPEN !!!\n"); // JMT
 581       continue;
 582     }
 583
 584     // -- Get sign of particle
 585     Float_t sign      = (particle->GetPdgCode() < 0) ? -1. : 1.;
 586
 587     // -- Get if particle is findable
 588     Float_t findable  = Float_t(fHelper->IsParticleFindable(idxMC));
 589
 590     // -- Get recStatus and pidStatus
 591     Float_t recStatus = 0.;
 592     Float_t recPid    = 0.;
 593
 594     // -- Get Reconstructed values
 595     Float_t etaRec = 0.;
 596     Float_t phiRec = 0.;
 597     Float_t ptRec  = 0.;
 598
 599     // -- Loop over all labels
 600     for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {
 601       if (idxMC == fLabelsRec[0][idxRec]) {
 602         recStatus = 1.;
 603
 604         if (idxMC == fLabelsRec[1][idxRec])
 605           recPid = 1.;
 606
 607         AliESDtrack *track = fESD->GetTrack(idxRec);
 608         if (track) {
 609           // if no track present (which should not happen)
 610           // -> pt = 0. , which is not in the looked at range
 611
 612           // -- Get Reconstructed eta/phi/pt
 613           etaRec = track->Eta();
 614           phiRec = track->Phi();
 615           ptRec  = track->Pt();
 616         }
 617         break;
 618       }
 619     } // for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {
 620
 621     // -- Fill THnSparse
 622     Double_t hnEff[15] = {fCentralityBin, particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), sign,
 623                           findable, recStatus, recPid, etaRec, phiRec, ptRec,
 624                           particle->Eta()-etaRec, particle->Phi()-phiRec, particle->Pt()-ptRec};
 625     fHnEff->Fill(hnEff);
 626
 627   } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
 628
 629   return;
 630 }
 631
 632 //________________________________________________________________________
 633 void AliAnalysisNetParticleEffCont::FillMCEffHistAOD() {
 634   // Fill efficiency THnSparse for AODs
 635
 636   Int_t nPart  = fArrayMC->GetEntriesFast();
 637
 638   for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
 639
 640     AliAODMCParticle* particle = static_cast<AliAODMCParticle*>(fArrayMC->At(idxMC));
 641
 642     // -- Check basic MC properties -> charged physical primary
 643     if (!fHelper->IsParticleAcceptedBasicCharged(particle))
 644       continue;
 645
 646     // -- Check rapidity window
 647     Double_t yMC;
 648     if (!fHelper->IsParticleAcceptedRapidity((TParticle*)particle, yMC))  // MW : das geht? wenn ja use c++ cast
 649       continue;
 650
 651     // -- Check if probeParticle / anti-probeParticle
 652     //    > skip check if PID is not required
 653     if (fHelper->GetUsePID() && TMath::Abs(particle->GetPdgCode()) != fPdgCode)
 654       continue;
 655
 656     // -- Get sign of particle
 657     Float_t sign      = (particle->GetPdgCode() < 0) ? -1. : 1.;
 658
 659     // -- Get if particle is findable
 660     Float_t findable  = 1.;// Float_t(fHelper->IsParticleFindable(idxMC)); XXX
 661
 662     // -- Get recStatus and pidStatus
 663     Float_t recStatus = 0.;
 664     Float_t recPid    = 0.;
 665
 666     // -- Get Reconstructed values
 667     Float_t etaRec = 0.;
 668     Float_t phiRec = 0.;
 669     Float_t ptRec  = 0.;
 670
 671     // -- Loop over all labels
 672     for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {
 673       if (idxMC == fLabelsRec[0][idxRec]) {
 674         recStatus = 1.;
 675
 676         if (idxMC == fLabelsRec[1][idxRec])
 677           recPid = 1.;
 678
 679         AliVTrack *track = NULL;
 680         if(fESD)
 681           track = fESD->GetTrack(idxRec);
 682         else if(fAOD)
 683           track = fAOD->GetTrack(idxRec);
 684
 685         if (track) {
 686           // if no track present (which should not happen)
 687           // -> pt = 0. , which is not in the looked at range
 688
 689           // -- Get Reconstructed eta/phi/pt
 690           etaRec = track->Eta();
 691           phiRec = track->Phi();
 692           ptRec  = track->Pt();
 693         }
 694         break;
 695       }
 696     } // for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {
 697
 698     // -- Fill THnSparse
 699     Double_t hnEff[15] = {fCentralityBin, particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), sign,
 700                           findable, recStatus, recPid, etaRec, phiRec, ptRec,
 701                           particle->Eta()-etaRec, particle->Phi()-phiRec, particle->Pt()-ptRec};
 702     fHnEff->Fill(hnEff);
 703
 704   } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
 705
 706   return;
 707 }