PWGCF/EBYE/NetParticle/AliAnalysisNetParticleDistribution.cxx

   1 //-*- Mode: C++ -*-
   2
   3 #include "TMath.h"
   4 #include "TAxis.h"
   5 #include "TSystem.h"
   6 #include "TProfile.h"
   7 #include "TH2F.h"
   8 #include "TH3F.h"
   9 #include "TFile.h"
  10 #include "TPRegexp.h"
  11
  12 #include "AliStack.h"
  13 #include "AliMCEvent.h"
  14 #include "AliMCParticle.h"
  15 #include "AliESDtrackCuts.h"
  16 #include "AliESDInputHandler.h"
  17 #include "AliESDpid.h"
  18 #include "AliCentrality.h"
  19 #include "AliTracker.h"
  20 #include "AliAODInputHandler.h"
  21 #include "AliAODEvent.h"
  22 #include "AliAODTrack.h"
  23
  24 #include "AliAnalysisNetParticleDistribution.h"
  25
  26 using namespace std;
  27
  28 // Task for NetParticle checks
  29 // Author: Jochen Thaeder <jochen@thaeder.de>
  30
  31 ClassImp(AliAnalysisNetParticleDistribution)
  32
  33 /*
  34  * ---------------------------------------------------------------------------------
  35  *                            Constructor / Destructor
  36  * ---------------------------------------------------------------------------------
  37  */
  38
  39 //________________________________________________________________________
  40 AliAnalysisNetParticleDistribution::AliAnalysisNetParticleDistribution() :
  41   fHelper(NULL),
  42
  43   fOutList(NULL),
  44
  45   fESDHandler(NULL),
  46   fPIDResponse(NULL),
  47   fESD(NULL),
  48   fAODHandler(NULL),
  49   fAOD(NULL),
  50   fIsMC(kFALSE),
  51   fMCEvent(NULL),
  52   fStack(NULL),
  53   fESDTrackCuts(NULL),
  54   fEtaMax(0.9),
  55   fPtRange(NULL),
  56   fAODtrackCutBit(1024),
  57   fNp(NULL),
  58   fNCorrNp(2),
  59   fCorrNp(NULL),
  60   fNMCNp(5),
  61   fMCNp(NULL),
  62   fNControlMCNp(5),
  63   fControlMCNp(NULL),
  64
  65   fHnTrackUnCorr(NULL) {
  66   // Constructor
  67
  68   AliLog::SetClassDebugLevel("AliAnalysisNetParticleDistribution",10);
  69 }
  70
  71 //________________________________________________________________________
  72 AliAnalysisNetParticleDistribution::~AliAnalysisNetParticleDistribution() {
  73   // Destructor
  74
  75   if (fNp) delete[] fNp;
  76
  77   for (Int_t ii = 0; ii < fNCorrNp; ++ii)
  78     if (fCorrNp[ii]) delete[] fCorrNp[ii];
  79   if (fCorrNp) delete[] fCorrNp;
  80
  81   for (Int_t ii = 0; ii < fNMCNp; ++ii)
  82     if (fMCNp[ii]) delete[] fMCNp[ii];
  83   if (fMCNp) delete[] fMCNp;
  84
  85   for (Int_t ii = 0; ii < fNControlMCNp; ++ii)
  86     if (fControlMCNp[ii]) delete[] fControlMCNp[ii];
  87   if (fControlMCNp) delete[] fControlMCNp;
  88
  89   return;
  90 }
  91
  92 /*
  93  * ---------------------------------------------------------------------------------
  94  *                                 Public Methods
  95  * ---------------------------------------------------------------------------------
  96  */
  97
  98 //________________________________________________________________________
  99 Int_t AliAnalysisNetParticleDistribution::Initialize(AliAnalysisNetParticleHelper* helper, AliESDtrackCuts* cuts, Bool_t isMC, Float_t *ptRange, Float_t etaMax, Int_t trackCutBit) {
 100   // -- Initialize
 101
 102   fHelper = helper;
 103   fESDTrackCuts = cuts;
 104   fIsMC = isMC;
 105   fPtRange = ptRange;
 106   fEtaMax = etaMax;
 107   fAODtrackCutBit = trackCutBit;
 108
 109   // ------------------------------------------------------------------
 110   // -- N particles / N anti-particles
 111   // ------------------------------------------------------------------
 112   //  Np          : arr[particle]
 113   //  CorrNp      : arr[corrSet][particle]
 114   //  MCNp        : arr[corrSet][particle] - MC
 115   //  ControlMCNp : arr[corrSet][particle] - Control MC
 116
 117   fNp     = new Float_t[2];
 118
 119   fCorrNp = new Float_t*[fNCorrNp];
 120   for (Int_t ii = 0 ; ii < fNCorrNp; ++ii)
 121     fCorrNp[ii] = new Float_t[2];
 122
 123   fMCNp = new Float_t*[fNMCNp];
 124   for (Int_t ii = 0 ; ii < fNMCNp; ++ii)
 125     fMCNp[ii] = new Float_t[2];
 126
 127   fControlMCNp = new Float_t*[fNControlMCNp];
 128   for (Int_t ii = 0 ; ii < fNControlMCNp; ++ii)
 129     fControlMCNp[ii] = new Float_t[2];
 130
 131   ResetEvent();
 132
 133   return 0;
 134 }
 135
 136 //________________________________________________________________________
 137 void AliAnalysisNetParticleDistribution::CreateHistograms(TList* outList) {
 138   // -- Add histograms to outlist
 139
 140   fOutList = outList;
 141
 142   // ------------------------------------------------------------------
 143   // -- Create net particle histograms
 144   // ------------------------------------------------------------------
 145
 146   AddHistSet("fHDist", "Uncorrected");
 147   AddHistSet("fHDistCorr0", "Corrected [without cross section correction]");
 148   AddHistSet("fHDistCorr1", "Corrected [with cross section correction]");
 149
 150   if (fIsMC) {
 151     AddHistSet("fHMCrapidity", "MC primary in |y| < 0.5");
 152     AddHistSet("fHMCptMin",    "MC primary in |y| + #it{p}_{T} > 0.1");
 153     AddHistSet("fHMCpt",       Form("MC primary in |y| < 0.5 + #it{p}_{T} [%.1f,%.1f]", fPtRange[0], fPtRange[1]));
 154     AddHistSet("fHMCeta",      Form("MC primary in |y| < 0.5 + |#eta| < %.1f", fEtaMax));
 155     AddHistSet("fHMCetapt",    Form("MC primary in |y| < 0.5 + |#eta| < %.1f + #it{p}_{T} [%.1f,%.1f]", fEtaMax, fPtRange[0], fPtRange[1]));
 156
 157     AddHistSet("fHControlMCLambdarapidity", "Control MC Lambda primary in |y| < 0.5");
 158     AddHistSet("fHControlMCLambdaptMin",    "Control MC Lambda primary in |y| + #it{p}_{T} > 0.1");
 159     AddHistSet("fHControlMCLambdapt",       Form("Control MC primary in |y| < 0.5 + #it{p}_{T} [%.1f,%.1f]", fPtRange[0], fPtRange[1]));
 160     AddHistSet("fHControlMCLambdaeta",      Form("Control MC primary in |y| < 0.5 + |#eta| < %.1f", fEtaMax));
 161     AddHistSet("fHControlMCLambdaetapt",    Form("Control MC primary in |y| < 0.5 + |#eta| < %.1f + #it{p}_{T} [%.1f,%.1f]", fEtaMax, fPtRange[0], fPtRange[1]));
 162   }
 163
 164   // ------------------------------------------------------------------
 165   // -- Get Probe Particle Container
 166   // ------------------------------------------------------------------
 167
 168   Double_t dCent[2] = {-0.5, 8.5};
 169   Int_t iCent       = 9;
 170
 171   Double_t dEta[2]  = {-0.9, 0.9}; // -> 37 bins
 172   Int_t iEta        = Int_t((dEta[1]-dEta[0]) / 0.075) +1 ;
 173
 174   Double_t dRap[2]  = {-0.5, 0.5};
 175   Int_t iRap        = Int_t((dRap[1]-dRap[0]) / 0.075) +1 ;
 176
 177   Double_t dPhi[2]  = {0.0, TMath::TwoPi()};
 178   Int_t iPhi        = 42;
 179
 180   Double_t dPt[2]   = {0.1, 3.0};
 181   Int_t iPt         = Int_t((dPt[1]-dPt[0]) / 0.075);
 182
 183   Double_t dSign[2] = {-1.5, 1.5};
 184   Int_t iSign       = 3;
 185
 186   //                          cent:     pt:     sign:     eta:     phi:       y
 187   Int_t    binShort[6] = {   iCent,    iPt,    iSign,    iEta,    iPhi,    iRap};
 188   Double_t minShort[6] = {dCent[0], dPt[0], dSign[0], dEta[0], dPhi[0], dRap[0]};
 189   Double_t maxShort[6] = {dCent[1], dPt[1], dSign[1], dEta[1], dPhi[1], dRap[1]};
 190
 191   // -- UnCorrected
 192   fOutList->Add(new THnSparseF("fHnTrackUnCorr", "cent:pt:sign:eta:phi:y", 6, binShort, minShort, maxShort));
 193   fHnTrackUnCorr = static_cast<THnSparseF*>(fOutList->Last());
 194   fHnTrackUnCorr->GetAxis(0)->SetTitle("centrality");
 195   fHnTrackUnCorr->GetAxis(1)->SetTitle("#it{p}_{T} (GeV/#it{c})");
 196   fHnTrackUnCorr->GetAxis(2)->SetTitle("sign");
 197   fHnTrackUnCorr->GetAxis(3)->SetTitle("#eta");
 198   fHnTrackUnCorr->GetAxis(4)->SetTitle("#varphi");
 199   fHnTrackUnCorr->GetAxis(5)->SetTitle("#it{y}");
 200
 201   fHelper->BinLogAxis(fHnTrackUnCorr, 1);
 202
 203   // ------------------------------------------------------------------
 204
 205   return;
 206 }
 207
 208 //________________________________________________________________________
 209 Int_t AliAnalysisNetParticleDistribution::SetupEvent(AliESDInputHandler *esdHandler, AliAODInputHandler *aodHandler,  AliMCEvent *mcEvent) {
 210   // -- Setup Event
 211
 212   ResetEvent();
 213
 214   // -- Get ESD objects
 215   if(esdHandler){
 216     fESDHandler  = esdHandler;
 217     fPIDResponse = esdHandler->GetPIDResponse();
 218     fESD         = fESDHandler->GetEvent();
 219   }
 220
 221   // -- Get AOD objects
 222   else if(aodHandler){
 223     fAODHandler  = aodHandler;
 224     fPIDResponse = aodHandler->GetPIDResponse();
 225     fAOD         = fAODHandler->GetEvent();
 226   }
 227
 228   // -- Get MC objects
 229   fMCEvent     = mcEvent;
 230   if (fMCEvent)
 231     fStack     = fMCEvent->Stack();
 232
 233   return 0;
 234 }
 235
 236 //________________________________________________________________________
 237 void AliAnalysisNetParticleDistribution::ResetEvent() {
 238   // -- Reset event
 239
 240   // -- Reset ESD Event
 241   fESD       = NULL;
 242
 243   // -- Reset AOD Event
 244   fAOD       = NULL;
 245
 246   // -- Reset MC Event
 247   if (fIsMC)
 248     fMCEvent = NULL;
 249
 250   // -- Reset N particles/anti-particles
 251   for (Int_t jj = 0; jj < 2; ++jj)
 252     fNp[jj] = 0.;
 253
 254   // -- Reset N corrected particles/anti-particles
 255   for (Int_t ii = 0; ii < fNCorrNp; ++ii)
 256     for (Int_t jj = 0; jj < 2; ++jj)
 257       fCorrNp[ii][jj] = 0.;
 258
 259   // -- Reset N MC particles/anti-particles
 260   for (Int_t ii = 0; ii < fNMCNp; ++ii)
 261     for (Int_t jj = 0; jj < 2; ++jj)
 262       fMCNp[ii][jj] = 0.;
 263
 264   // -- Reset N control MC particles/anti-particles
 265   for (Int_t ii = 0; ii < fNControlMCNp; ++ii)
 266     for (Int_t jj = 0; jj < 2; ++jj)
 267       fControlMCNp[ii][jj] = 0.;
 268 }
 269
 270 //________________________________________________________________________
 271 Int_t AliAnalysisNetParticleDistribution::Process() {
 272   // -- Process NetParticle Distributions
 273
 274   // -- Fill ESD tracks
 275   if (fESD) ProcessESDTracks();
 276
 277   // -- Fill AOD tracks
 278   else if (fAOD) ProcessAODTracks();
 279
 280   // -- Fill MC truth particles (missing for AOD XXX)
 281   if (fIsMC)  {
 282     ProcessStackParticles();
 283     ProcessStackControlParticles();
 284   }
 285
 286   return 0;
 287 }
 288
 289 //________________________________________________________________________
 290 Int_t AliAnalysisNetParticleDistribution::ProcessESDTracks() {
 291   // -- Process ESD tracks and fill histograms
 292
 293   for (Int_t idxTrack = 0; idxTrack < fESD->GetNumberOfTracks(); ++idxTrack) {
 294     AliESDtrack *track = fESD->GetTrack(idxTrack);
 295
 296     // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 297     // -- Check track cuts
 298     // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 299
 300     // -- Check if charged track is accepted for basic parameters
 301     if (!fHelper->IsTrackAcceptedBasicCharged(track))
 302       continue;
 303
 304     // -- Check if accepted
 305     if (!fESDTrackCuts->AcceptTrack(track))
 306       continue;
 307
 308     // -- Check if accepted in rapidity window
 309     Double_t yP;
 310     if (!fHelper->IsTrackAcceptedRapidity(track, yP))
 311       continue;
 312
 313     // -- Check if accepted bt PID from TPC or TPC+TOF
 314     Double_t pid[2];
 315     if (!fHelper->IsTrackAcceptedPID(track, pid))
 316       continue;
 317
 318     // -- Check if accepted with thighter DCA cuts
 319     if (fHelper->IsTrackAcceptedDCA(track))
 320       continue;
 321
 322     // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 323     // -- Fill Probe Particle
 324     // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 325
 326     Double_t aTrack[6] = {
 327       Double_t(fHelper->GetCentralityBin()),       //  0 centrality
 328       track->Pt(),                    //  1 pt
 329       track->GetSign(),               //  2 sign
 330       track->Eta(),                   //  3 eta
 331       track->Phi(),                   //  4 phi
 332       yP                              //  5 rapidity
 333     };
 334
 335     fHnTrackUnCorr->Fill(aTrack);
 336
 337     // -- Count particle / anti-particle
 338     // ------------------------------------------------------------------
 339     //  idxPart = 0 -> anti particle
 340     //  idxPart = 1 -> particle
 341
 342     Int_t idxPart = (track->GetSign() < 0) ? 0 : 1;
 343     fNp[idxPart] += 1.;
 344
 345     for (Int_t ii = 0; ii < fNCorrNp; ++ii)
 346       fCorrNp[ii][idxPart] += fHelper->GetTrackbyTrackCorrectionFactor(aTrack, ii);
 347
 348   } // for (Int_t idxTrack = 0; idxTrack < fESD->GetNumberOfTracks(); ++idxTrack) {
 349
 350   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 351   // -- Fill Particle Fluctuation Histograms
 352   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 353
 354   // -- Uncorrected
 355   FillHistSet("fHDist", fNp);
 356
 357   // -- Corrected
 358   for (Int_t ii = 0 ; ii < fNCorrNp; ++ii)
 359     FillHistSet(Form("fHDistCorr%d", ii), fCorrNp[ii]);
 360
 361   return 0;
 362 }
 363
 364 //________________________________________________________________________
 365 Int_t AliAnalysisNetParticleDistribution::ProcessAODTracks() {
 366   // -- Process ESD tracks and fill histograms
 367
 368   for (Int_t idxTrack = 0; idxTrack < fAOD->GetNumberOfTracks(); ++idxTrack) {
 369     AliAODTrack *track = (AliAODTrack*)fAOD->GetTrack(idxTrack);
 370
 371     // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 372     // -- Check track cuts
 373     // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 374
 375     // -- Check if charged track is accepted for basic parameters
 376     if (!fHelper->IsTrackAcceptedBasicCharged(track))
 377       continue;
 378
 379     // -- Check if accepted
 380     if(!track->TestFilterBit(fAODtrackCutBit))
 381       continue;
 382
 383     // -- Check if accepted in rapidity window
 384     Double_t yP;
 385     if (!fHelper->IsTrackAcceptedRapidity(track, yP))
 386       continue;
 387
 388     // -- Check if accepted bt PID from TPC or TPC+TOF
 389     Double_t pid[2];
 390     if (!fHelper->IsTrackAcceptedPID(track, pid))
 391       continue;
 392
 393     // -- Check if accepted with thighter DCA cuts XXX
 394     // if (fHelper->IsTrackAcceptedDCA(track))
 395     //  continue;
 396
 397     // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 398     // -- Fill Probe Particle
 399     // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 400
 401     Double_t aTrack[6] = {
 402       Double_t(fHelper->GetCentralityBin()),       //  0 centrality
 403       track->Pt(),                    //  1 pt
 404       track->Charge(),                //  2 sign
 405       track->Eta(),                   //  3 eta
 406       track->Phi(),                   //  4 phi
 407       yP                              //  5 rapidity
 408     };
 409
 410     fHnTrackUnCorr->Fill(aTrack);
 411
 412     // -- Count particle / anti-particle
 413     // ------------------------------------------------------------------
 414     //  idxPart = 0 -> anti particle
 415     //  idxPart = 1 -> particle
 416
 417     Int_t idxPart = (track->Charge() < 0) ? 0 : 1;
 418     fNp[idxPart] += 1.;
 419
 420     for (Int_t ii = 0; ii < fNCorrNp; ++ii)
 421       fCorrNp[ii][idxPart] += fHelper->GetTrackbyTrackCorrectionFactor(aTrack, ii);
 422
 423   } // for (Int_t idxTrack = 0; idxTrack < fESD->GetNumberOfTracks(); ++idxTrack) {
 424
 425   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 426   // -- Fill Particle Fluctuation Histograms
 427   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 428
 429   // -- Uncorrected
 430   FillHistSet("fHDist", fNp);
 431
 432   // -- Corrected
 433   for (Int_t ii = 0 ; ii < fNCorrNp; ++ii)
 434     FillHistSet(Form("fHDistCorr%d", ii), fCorrNp[ii]);
 435
 436   return 0;
 437 }
 438
 439 //________________________________________________________________________
 440 Int_t AliAnalysisNetParticleDistribution::ProcessStackParticles() {
 441   // -- Process primary particles from the stack and fill histograms
 442
 443   Int_t pdgCode    = AliPID::ParticleCode(fHelper->GetParticleSpecies());
 444
 445   for (Int_t idxMC = 0; idxMC < fStack->GetNprimary(); ++idxMC) {
 446     TParticle* particle = fStack->Particle(idxMC);
 447     if (!particle)
 448       continue;
 449
 450     // -- Check basic MC properties -> charged physical primary
 451     if (!fHelper->IsParticleAcceptedBasicCharged(particle, idxMC))
 452       continue;
 453
 454     // -- Check if particle / anti-particle
 455     if (TMath::Abs(particle->GetPdgCode()) != pdgCode)
 456       continue;
 457
 458     // -- Get particle : 0 anti-particle / 1 particle
 459     Int_t idxPart = (particle->GetPdgCode() < 0) ? 0 : 1;
 460
 461     // >> NOW only anti-particle / particle
 462     // >> With idxPart
 463
 464     // -- Check rapidity window
 465     Double_t yMC;
 466     if (!fHelper->IsParticleAcceptedRapidity(particle, yMC))
 467       continue;
 468
 469     fMCNp[0][idxPart] += 1.;         // -> MCrapidity
 470
 471     // -- Check acceptance
 472     if (particle->Pt() > 0.1 )
 473       fMCNp[1][idxPart] += 1.;       // -> MCptMin
 474
 475     if (particle->Pt() > fPtRange[0] && particle->Pt() <= fPtRange[1])
 476       fMCNp[2][idxPart] += 1.;       // -> MCpt
 477
 478     if (TMath::Abs(particle->Eta()) <= fEtaMax)
 479       fMCNp[3][idxPart] += 1.;       // -> MCeta
 480
 481     if (particle->Pt() > fPtRange[0] && particle->Pt() > fPtRange[1] && TMath::Abs(particle->Eta()) <= fEtaMax)
 482       fMCNp[4][idxPart] += 1.;       // -> MCetapt
 483
 484   } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
 485
 486   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 487   // -- Fill Particle Fluctuation Histograms
 488   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 489
 490   FillHistSet("fHMCrapidity", fMCNp[0]);
 491   FillHistSet("fHMCptMin",    fMCNp[1]);
 492   FillHistSet("fHMCpt",       fMCNp[2]);
 493   FillHistSet("fHMCeta",      fMCNp[3]);
 494   FillHistSet("fHMCetapt",    fMCNp[4]);
 495
 496   return 0;
 497 }
 498
 499 //________________________________________________________________________
 500 Int_t AliAnalysisNetParticleDistribution::ProcessStackControlParticles() {
 501   // -- Process primary particles from the stack and fill histograms
 502
 503   Int_t pdgCode    = fHelper->GetControlParticleCode();
 504   Bool_t isNeutral  = fHelper->IsControlParticleNeutral();
 505   const Char_t* name = fHelper->GetControlParticleName().Data();
 506
 507   for (Int_t idxMC = 0; idxMC < fStack->GetNprimary(); ++idxMC) {
 508     TParticle* particle = fStack->Particle(idxMC);
 509     if (!particle)
 510       continue;
 511
 512     // -- Check basic MC properties -> neutral or charged physical primary
 513     if (isNeutral) {
 514       if (!fHelper->IsParticleAcceptedBasicNeutral(particle, idxMC))
 515         continue;
 516     }
 517     else {
 518       if (!fHelper->IsParticleAcceptedBasicCharged(particle, idxMC))
 519         continue;
 520     }
 521
 522     // -- Check if particle / anti-particle
 523     if (TMath::Abs(particle->GetPdgCode()) != pdgCode)
 524       continue;
 525
 526     // -- Get particle : 0 anti-particle / 1 particle
 527     Int_t idxPart = (particle->GetPdgCode() < 0) ? 0 : 1;
 528
 529     // >> NOW only anti-particle / particle
 530     // >> With idxPart
 531
 532     // -- Check rapidity window
 533     Double_t yMC;
 534     if (!fHelper->IsParticleAcceptedRapidity(particle, yMC))
 535       continue;
 536
 537     fControlMCNp[0][idxPart] += 1.;         // -> ControlMCrapidity
 538
 539     // -- Check acceptance
 540     if (particle->Pt() > 0.1 )
 541       fControlMCNp[1][idxPart] += 1.;       // -> ControlMCptMin
 542
 543     if (particle->Pt() > fPtRange[0] && particle->Pt() <= fPtRange[1])
 544       fControlMCNp[2][idxPart] += 1.;       // -> ControlMCpt
 545
 546     if (TMath::Abs(particle->Eta()) <= fEtaMax)
 547       fControlMCNp[3][idxPart] += 1.;       // -> ControlMCeta
 548
 549     if (particle->Pt() > fPtRange[0] && particle->Pt() > fPtRange[1] && TMath::Abs(particle->Eta()) <= fEtaMax)
 550       fControlMCNp[4][idxPart] += 1.;       // -> ControlMCetapt
 551   } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
 552
 553   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 554   // -- Fill Particle Fluctuation Histograms
 555   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 556
 557   FillHistSet(Form("fHControlMC%srapidity", name), fControlMCNp[0], 0);
 558   FillHistSet(Form("fHControlMC%sptMin",    name), fControlMCNp[1], 1);
 559   FillHistSet(Form("fHControlMC%spt",       name), fControlMCNp[2], 2);
 560   FillHistSet(Form("fHControlMC%seta",      name), fControlMCNp[3], 3);
 561   FillHistSet(Form("fHControlMC%setapt",    name), fControlMCNp[4], 4);
 562
 563   return 0;
 564 }
 565
 566 ////////////////////////////////////////////////////////////////////////////////////////////
 567
 568
 569 //________________________________________________________________________
 570 void  AliAnalysisNetParticleDistribution::AddHistSet(const Char_t *name, const Char_t *title)  {
 571   // -- Add histogram sets for particle and anti-particle
 572
 573   const Char_t* aPartNames[]             = {"pbar",          "p"};
 574   const Char_t* aPartTitles[]            = {"Anti-Proton",   "Proton"};
 575   const Char_t* aPartAxisTitles[]        = {"#bar{p}",       "p"};
 576
 577   const Char_t* aControlPartNames[]      = {"lambdabar",     "lambda"};
 578   const Char_t* aControlPartTitles[]     = {"Anti-Lambda",   "Lambda"};
 579
 580   fOutList->Add(new TList);
 581   TList *list = static_cast<TList*>(fOutList->Last());
 582   list->SetName(name) ;
 583   list->SetOwner(kTRUE);
 584
 585   TString sName(name);
 586
 587   TString sPtTitle("");
 588   if (!sName.Contains("fHMC") || !sName.Contains("fHControlMC"))
 589     sPtTitle += Form("#it{p}_{T} [%.1f,%.1f]", fPtRange[0], fPtRange[1]);
 590
 591   for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
 592     list->Add(new TH2F(Form("%s%s", name, aPartNames[idxPart]),
 593                        Form("%s %s Dist %s;Centrality;N_{%s}", title, aPartTitles[idxPart], sPtTitle.Data(), aPartAxisTitles[idxPart]),
 594                        24, -0.5, 11.5, 501, -0.5, 500.49));
 595   } // for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
 596
 597   list->Add(new TH2F(Form("%sNet%s", name, aPartNames[1]),
 598                      Form("%s Net %s Dist %s;Centrality;N_{p} - N_{#bar{p}}", title, aPartTitles[1], sPtTitle.Data()), 24, -0.5, 11.5, 201, -100.5, 100.49));
 599
 600   list->Add(new TProfile(Form("%sNet%sM", name, aPartNames[1]),
 601                          Form("%s Net %s Dist %s;Centrality;N_{p} - N_{#bar{p}}", title, aPartTitles[1], sPtTitle.Data()),  24, -0.5, 11.5));
 602   list->Add(new TProfile(Form("%sNet%s2M", name, aPartNames[1]),
 603                          Form("%s (Net %s)^{2} Dist %s;Centrality;(N_{p} - N_{#bar{p}})^{2}", title, aPartTitles[1], sPtTitle.Data()), 24, -0.5, 11.5));
 604   list->Add(new TProfile(Form("%sNet%s3M", name, aPartNames[1]),
 605                          Form("%s (Net %s)^{3} Dist %s;Centrality;(N_{p} - N_{#bar{p}})^{3}", title, aPartTitles[1], sPtTitle.Data()), 24, -0.5, 11.5));
 606   list->Add(new TProfile(Form("%sNet%s4M", name, aPartNames[1]),
 607                          Form("%s (Net %s)^{4} Dist %s;Centrality;(N_{p} - N_{#bar{p}})^{4}", title, aPartTitles[1], sPtTitle.Data()), 24, -0.5, 11.5));
 608   list->Add(new TProfile(Form("%sNet%s5M", name, aPartNames[1]),
 609                          Form("%s (Net %s)^{5} Dist %s;Centrality;(N_{p} - N_{#bar{p}})^{5}", title, aPartTitles[1], sPtTitle.Data()), 24, -0.5, 11.5));
 610   list->Add(new TProfile(Form("%sNet%s6M", name, aPartNames[1]),
 611                          Form("%s (Net %s)^{6} Dist %s;Centrality;(N_{p} - N_{#bar{p}})^{6}", title, aPartTitles[1], sPtTitle.Data()), 24, -0.5, 11.5));
 612
 613   list->Add(new TH2F(Form("%sNet%sOverSum", name, aPartNames[1]),
 614                      Form("%s (Net %s)/ Sum Dist %s;Centrality;(N_{p} - N_{#bar{p}})/(N_{p} + N_{#bar{p}})", title, aPartTitles[1], sPtTitle.Data()),
 615                      24, -0.5, 11.5, 801, -50.5, 50.49));
 616
 617   if (sName.Contains("fHControlMC")) {
 618     for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
 619       list->Add(new TH2F(Form("%s%sOver%s", name, aPartNames[idxPart], aControlPartNames[idxPart]),
 620                          Form("%s %s / %s Dist %s;Centrality;N_{%s}/N_{Tracks}",
 621                               title, aPartTitles[idxPart], aControlPartTitles[idxPart], sPtTitle.Data(), aPartAxisTitles[idxPart]),
 622                          24, -0.5, 11.5, 101, 0., 1.));
 623     } // for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
 624   }
 625
 626   //  if (sName.Contains("fHDist")) {
 627   //    list->Add(new TH3F(Form("%sNet%s_RecVsMC", name, aPartNames[1]),
 628   //                   Form("%s Net %s Dist - Rec Vs MC - %s;Centrality;(N_{p} - N_{#bar{p}})_{rec};(N_{p} - N_{#bar{p}})_{MC}",
 629   //                        title, aPartTitles[1], sPtTitle.Data()),
 630   //                   24, -0.5, 11.5, 201, -100.5, 100.49, 201, -100.5, 100.49));
 631   //  }
 632
 633   return;
 634 }
 635
 636 //________________________________________________________________________
 637 void AliAnalysisNetParticleDistribution::FillHistSet(const Char_t *name, Float_t *np, Int_t controlIdx)  {
 638   // -- Add histogram sets for particle and anti-particle
 639
 640   TList *list = static_cast<TList*>(fOutList->FindObject(name));
 641
 642   Float_t centralityBin = fHelper->GetCentralityBin();
 643
 644   (static_cast<TH2F*>(list->FindObject(Form("%spbar", name))))->Fill(centralityBin, np[0]);
 645   (static_cast<TH2F*>(list->FindObject(Form("%sp",    name))))->Fill(centralityBin, np[1]);
 646
 647   Float_t sumNp    = np[1]+np[0];
 648   Float_t deltaNp  = np[1]-np[0];
 649   Float_t deltaNp2 = deltaNp * deltaNp;
 650   Float_t deltaNp3 = deltaNp2 * deltaNp;
 651
 652   (static_cast<TH2F*>(list->FindObject(Form("%sNetp",  name))))->Fill(centralityBin, deltaNp);
 653
 654   (static_cast<TProfile*>(list->FindObject(Form("%sNetpM",  name))))->Fill(centralityBin, deltaNp);
 655   (static_cast<TProfile*>(list->FindObject(Form("%sNetp2M", name))))->Fill(centralityBin, deltaNp2);
 656   (static_cast<TProfile*>(list->FindObject(Form("%sNetp3M", name))))->Fill(centralityBin, deltaNp2*deltaNp);
 657   (static_cast<TProfile*>(list->FindObject(Form("%sNetp4M", name))))->Fill(centralityBin, deltaNp2*deltaNp2);
 658   (static_cast<TProfile*>(list->FindObject(Form("%sNetp5M", name))))->Fill(centralityBin, deltaNp3*deltaNp2);
 659   (static_cast<TProfile*>(list->FindObject(Form("%sNetp6M", name))))->Fill(centralityBin, deltaNp3*deltaNp3);
 660
 661   (static_cast<TH2F*>(list->FindObject(Form("%sNetpOverSum", name))))->Fill(centralityBin, deltaNp/sumNp);
 662
 663   TString sName(name);
 664   if (sName.Contains("fHControlMC") && controlIdx >= 0) {
 665     (static_cast<TH2F*>(list->FindObject(Form("%spbarOverlambdabar", name))))->Fill(centralityBin, np[0]/fControlMCNp[controlIdx][0]);
 666     (static_cast<TH2F*>(list->FindObject(Form("%spOverlambda",       name))))->Fill(centralityBin, np[1]/fControlMCNp[controlIdx][1]);
 667   }
 668
 669   //  if (sName.Contains("fHDist")) {
 670   //    Float_t deltaMCNp  = fMCNp[controlIdx][1]-fMCNp[controldx][0];
 671   //    (static_cast<TH3F*>(list->FindObject(Form("%sNetp_RecVsMC", name))))->Fill(centralityBin, deltaNp, deltaMCNp);
 672   //  }
 673
 674   return;
 675 }
 676