PWGLF/SPECTRA/ChargedHadrons/dNdPt/AlidNdPtAnalysisPbPbAOD.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   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 // AlidNdPtAnalysisPbPbAOD class.
  17 //
  18 // Author: P. Luettig, 15.05.2013
  19 // last modified: 17.10.2013
  20 //------------------------------------------------------------------------------
  21
  22
  23 #include "AlidNdPtAnalysisPbPbAOD.h"
  24
  25
  26 using namespace std;
  27
  28 ClassImp(AlidNdPtAnalysisPbPbAOD)
  29
  30
  31
  32 AlidNdPtAnalysisPbPbAOD::AlidNdPtAnalysisPbPbAOD(const char *name) : AliAnalysisTaskSE(name),
  33 fOutputList(0),
  34 // Histograms
  35 fPt(0),
  36 fMCPt(0),
  37 fZvPtEtaCent(0),
  38 fPhiPtEtaCent(0),
  39 fMCRecPrimZvPtEtaCent(0),
  40 fMCGenZvPtEtaCent(0),
  41 fMCRecSecZvPtEtaCent(0),
  42 fMCRecPrimPhiPtEtaCent(0),
  43 fMCGenPhiPtEtaCent(0),
  44 fMCRecSecPhiPtEtaCent(0),
  45 fEventStatistics(0),
  46 fEventStatisticsCentrality(0),
  47 fMCEventStatisticsCentrality(0),
  48 fAllEventStatisticsCentrality(0),
  49 fEventStatisticsCentralityTrigger(0),
  50 fZvMultCent(0),
  51 fTriggerStatistics(0),
  52 fMCTrackPdgCode(0),
  53 fMCTrackStatusCode(0),
  54 fCharge(0),
  55 fMCCharge(0),
  56 fMCPdgPt(0),
  57 fMCHijingPrim(0),
  58 fDCAPtAll(0),
  59 fDCAPtAccepted(0),
  60 fMCDCAPtSecondary(0),
  61 fMCDCAPtPrimary(0),
  62 fCutPercClusters(0),
  63 fCutPercCrossed(0),
  64 fCrossCheckRowsLength(0),
  65 fCrossCheckClusterLength(0),
  66 fCrossCheckRowsLengthAcc(0),
  67 fCrossCheckClusterLengthAcc(0),
  68 //global
  69 fIsMonteCarlo(0),
  70 // event cut variables
  71 fCutMaxZVertex(10.),
  72 // track kinematic cut variables
  73 fCutPtMin(0.15),
  74 fCutPtMax(200.),
  75 fCutEtaMin(-0.8),
  76 fCutEtaMax(0.8),
  77 // track quality cut variables
  78 fUseRelativeCuts(kFALSE),
  79 fCutRequireTPCRefit(kTRUE),
  80 fCutMinNumberOfClusters(60),
  81 fCutPercMinNumberOfClusters(0.2),
  82 fCutMinNumberOfCrossedRows(120.),
  83 fCutPercMinNumberOfCrossedRows(0.2),
  84 fCutMinRatioCrossedRowsOverFindableClustersTPC(0.8),
  85 fCutMaxChi2PerClusterTPC(4.),
  86 fCutMaxFractionSharedTPCClusters(0.4),
  87 fCutMaxDCAToVertexZ(3.0),
  88 fCutMaxDCAToVertexXY(3.0),
  89 fCutRequireITSRefit(kTRUE),
  90 fCutMaxChi2PerClusterITS(36.),
  91 fCutDCAToVertex2D(kFALSE),
  92 fCutRequireSigmaToVertex(kFALSE),
  93 fCutMaxDCAToVertexXYPtDepPar0(0.0182),
  94 fCutMaxDCAToVertexXYPtDepPar1(0.0350),
  95 fCutMaxDCAToVertexXYPtDepPar2(1.01),
  96 fCutAcceptKinkDaughters(kFALSE),
  97 fCutMaxChi2TPCConstrainedGlobal(36.),
  98 fCutLengthInTPCPtDependent(kFALSE),
  99 fPrefactorLengthInTPCPtDependent(1),
 100 // binning for THnSparse
 101 fMultNbins(0),
 102 fPtNbins(0),
 103 fPtCorrNbins(0),
 104 fPtCheckNbins(0),
 105 fEtaNbins(0),
 106 fEtaCheckNbins(0),
 107 fZvNbins(0),
 108 fCentralityNbins(0),
 109 fPhiNbins(0),
 110 fBinsMult(0),
 111 fBinsPt(0),
 112 fBinsPtCorr(0),
 113 fBinsPtCheck(0),
 114 fBinsEta(0),
 115 fBinsEtaCheck(0),
 116 fBinsZv(0),
 117 fBinsCentrality(0),
 118 fBinsPhi(0)
 119 {
 120
 121   for(Int_t i = 0; i < cqMax; i++)
 122   {
 123     fCrossCheckAll[i] = 0;
 124     fCrossCheckAcc[i] = 0;
 125   }
 126
 127   fMultNbins = 0;
 128   fPtNbins = 0;
 129   fPtCorrNbins = 0;
 130   fPtCheckNbins = 0;
 131   fEtaNbins = 0;
 132   fEtaCheckNbins = 0;
 133   fZvNbins = 0;
 134   fCentralityNbins = 0;
 135   fBinsMult = 0;
 136   fBinsPt = 0;
 137   fBinsPtCorr = 0;
 138   fBinsPtCheck = 0;
 139   fBinsEta = 0;
 140   fBinsEtaCheck = 0;
 141   fBinsZv = 0;
 142   fBinsCentrality = 0;
 143   fBinsPhi = 0;
 144
 145   DefineOutput(1, TList::Class());
 146 }
 147
 148 // destructor
 149 AlidNdPtAnalysisPbPbAOD::~AlidNdPtAnalysisPbPbAOD()
 150 {
 151
 152   if(fZvPtEtaCent) delete fZvPtEtaCent; fZvPtEtaCent = 0;
 153   if(fPt) delete fPt; fPt = 0;
 154   if(fMCRecPrimZvPtEtaCent) delete fMCRecPrimZvPtEtaCent; fMCRecPrimZvPtEtaCent = 0;
 155   if(fMCGenZvPtEtaCent) delete fMCGenZvPtEtaCent; fMCGenZvPtEtaCent = 0;
 156   if(fMCRecSecZvPtEtaCent) delete fMCRecSecZvPtEtaCent; fMCRecSecZvPtEtaCent = 0;
 157   if(fMCPt) delete fMCPt; fMCPt = 0;
 158   if(fEventStatistics) delete fEventStatistics; fEventStatistics = 0;
 159   if(fEventStatisticsCentrality) delete fEventStatisticsCentrality; fEventStatisticsCentrality = 0;
 160   if(fMCEventStatisticsCentrality) delete fMCEventStatisticsCentrality; fMCEventStatisticsCentrality = 0;
 161   if(fAllEventStatisticsCentrality) delete fAllEventStatisticsCentrality; fAllEventStatisticsCentrality = 0;
 162   if(fEventStatisticsCentralityTrigger) delete fEventStatisticsCentralityTrigger; fEventStatisticsCentralityTrigger = 0;
 163   if(fZvMultCent) delete fZvMultCent; fZvMultCent = 0;
 164   if(fTriggerStatistics) delete fTriggerStatistics; fTriggerStatistics = 0;
 165   if(fMCTrackPdgCode) delete fMCTrackPdgCode; fMCTrackPdgCode = 0;
 166   if(fMCTrackStatusCode) delete fMCTrackStatusCode; fMCTrackStatusCode = 0;
 167   if(fCharge) delete fCharge; fCharge = 0;
 168   if(fMCCharge) delete fMCCharge; fMCCharge = 0;
 169   if(fMCPdgPt) delete fMCPdgPt; fMCPdgPt = 0;
 170   if(fMCHijingPrim) delete fMCHijingPrim; fMCHijingPrim = 0;
 171   if(fDCAPtAll) delete fDCAPtAll; fDCAPtAll = 0;
 172   if(fDCAPtAccepted) delete fDCAPtAccepted; fDCAPtAccepted = 0;
 173   if(fMCDCAPtSecondary) delete fMCDCAPtSecondary; fMCDCAPtSecondary = 0;
 174   if(fMCDCAPtPrimary) delete fMCDCAPtPrimary; fMCDCAPtPrimary = 0;
 175   if(fMCDCAPtSecondary) delete fMCDCAPtSecondary; fMCDCAPtSecondary = 0;
 176   if(fMCDCAPtPrimary) delete fMCDCAPtPrimary; fMCDCAPtPrimary = 0;
 177
 178   for(Int_t i = 0; i < cqMax; i++)
 179   {
 180     if(fCrossCheckAll[i]) delete fCrossCheckAll[i]; fCrossCheckAll[i] = 0;
 181     if(fCrossCheckAcc[i]) delete fCrossCheckAcc[i]; fCrossCheckAcc[i] = 0;
 182   }
 183
 184 }
 185
 186 void AlidNdPtAnalysisPbPbAOD::UserCreateOutputObjects()
 187 {
 188   // create all output histograms here
 189   OpenFile(1, "RECREATE");
 190
 191   fOutputList = new TList();
 192   fOutputList->SetOwner();
 193
 194   //define default binning
 195   Double_t binsMultDefault[48] = {-0.5, 0.5 , 1.5 , 2.5 , 3.5 , 4.5 , 5.5 , 6.5 , 7.5 , 8.5,9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5,19.5, 20.5, 30.5, 40.5 , 50.5 , 60.5 , 70.5 , 80.5 , 90.5 , 100.5,200.5, 300.5, 400.5, 500.5, 600.5, 700.5, 800.5, 900.5, 1000.5, 2000.5, 3000.5, 4000.5, 5000.5, 6000.5, 7000.5, 8000.5, 9000.5, 10000.5 };
 196   Double_t binsPtDefault[82] = {0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0, 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 110.0, 120.0, 130.0, 140.0, 150.0, 160.0, 180.0, 200.0};
 197   Double_t binsPtCorrDefault[37] = {0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 3.0, 4.0, 200.0};
 198   Double_t binsEtaDefault[31] = {-1.5,-1.4,-1.3,-1.2,-1.1,-1.0,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5};
 199   Double_t binsZvDefault[13] = {-30.,-25.,-20.,-15.,-10.,-5.,0.,5.,10.,15.,20.,25.,30.};
 200   Double_t binsCentralityDefault[12] = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.};
 201
 202   Double_t binsPhiDefault[37] = { 0, 0.174533, 0.349066, 0.523599, 0.698132, 0.872665, 1.0472, 1.22173, 1.39626, 1.5708, 1.74533, 1.91986, 2.0944, 2.26893, 2.44346, 2.61799, 2.79253, 2.96706, 3.14159, 3.31613, 3.49066, 3.66519, 3.83972, 4.01426, 4.18879, 4.36332, 4.53786, 4.71239, 4.88692, 5.06145, 5.23599, 5.41052, 5.58505, 5.75959, 5.93412, 6.10865, 2.*TMath::Pi()};
 203
 204   Double_t binsPtCheckDefault[20] = {0.,0.15,0.5,1.0,2.0,3.0,4.0, 5.0, 10.0, 13.0, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0, 70.0, 100.0, 150.0, 200.0};
 205   Double_t binsEtaCheckDefault[7] = {-1.0,-0.8,-0.4,0.,0.4,0.8,1.0};
 206
 207   // if no binning is set, use the default
 208   if (!fBinsMult)       { SetBinsMult(48,binsMultDefault); }
 209   if (!fBinsPt)         { SetBinsPt(82,binsPtDefault); }
 210   if (!fBinsPtCorr)     { SetBinsPtCorr(37,binsPtCorrDefault); }
 211   if (!fBinsPtCheck)    { SetBinsPtCheck(20,binsPtCheckDefault); }
 212   if (!fBinsEta)        { SetBinsEta(31,binsEtaDefault); }
 213   if (!fBinsEtaCheck)   { SetBinsEtaCheck(7,binsEtaCheckDefault); }
 214   if (!fBinsZv)         { SetBinsZv(13,binsZvDefault); }
 215   if (!fBinsCentrality) { SetBinsCentrality(12,binsCentralityDefault); }
 216   if (!fBinsPhi)        { SetBinsPhi(37,binsPhiDefault); }
 217
 218   Int_t binsZvPtEtaCent[4]={fZvNbins-1,fPtNbins-1,fEtaNbins-1,fCentralityNbins-1};
 219   Int_t binsPhiPtEtaCent[4]={fPhiNbins-1,fPtNbins-1,fEtaNbins-1,fCentralityNbins-1};
 220   Int_t binsZvMultCent[3]={fZvNbins-1,fMultNbins-1,fCentralityNbins-1};
 221
 222   // define Histograms
 223   fZvPtEtaCent = new THnSparseF("fZvPtEtaCent","Zv:Pt:Eta:Centrality",4,binsZvPtEtaCent);
 224   fZvPtEtaCent->SetBinEdges(0,fBinsZv);
 225   fZvPtEtaCent->SetBinEdges(1,fBinsPt);
 226   fZvPtEtaCent->SetBinEdges(2,fBinsEta);
 227   fZvPtEtaCent->SetBinEdges(3,fBinsCentrality);
 228   fZvPtEtaCent->GetAxis(0)->SetTitle("Zv (cm)");
 229   fZvPtEtaCent->GetAxis(1)->SetTitle("Pt (GeV/c)");
 230   fZvPtEtaCent->GetAxis(2)->SetTitle("Eta");
 231   fZvPtEtaCent->GetAxis(3)->SetTitle("Centrality");
 232   fZvPtEtaCent->Sumw2();
 233
 234   fPhiPtEtaCent = new THnSparseF("fPhiPtEtaCent","Phi:Pt:Eta:Centrality",4,binsPhiPtEtaCent);
 235   fPhiPtEtaCent->SetBinEdges(0,fBinsPhi);
 236   fPhiPtEtaCent->SetBinEdges(1,fBinsPt);
 237   fPhiPtEtaCent->SetBinEdges(2,fBinsEta);
 238   fPhiPtEtaCent->SetBinEdges(3,fBinsCentrality);
 239   fPhiPtEtaCent->GetAxis(0)->SetTitle("Phi (cm)");
 240   fPhiPtEtaCent->GetAxis(1)->SetTitle("Pt (GeV/c)");
 241   fPhiPtEtaCent->GetAxis(2)->SetTitle("Eta");
 242   fPhiPtEtaCent->GetAxis(3)->SetTitle("Centrality");
 243   fPhiPtEtaCent->Sumw2();
 244
 245
 246
 247   fMCRecPrimZvPtEtaCent = new THnSparseF("fMCRecPrimZvPtEtaCent","mcZv:mcPt:mcEta:Centrality",4,binsZvPtEtaCent);
 248   fMCRecPrimZvPtEtaCent->SetBinEdges(0,fBinsZv);
 249   fMCRecPrimZvPtEtaCent->SetBinEdges(1,fBinsPt);
 250   fMCRecPrimZvPtEtaCent->SetBinEdges(2,fBinsEta);
 251   fMCRecPrimZvPtEtaCent->SetBinEdges(3,fBinsCentrality);
 252   fMCRecPrimZvPtEtaCent->GetAxis(0)->SetTitle("MC Zv (cm)");
 253   fMCRecPrimZvPtEtaCent->GetAxis(1)->SetTitle("MC Pt (GeV/c)");
 254   fMCRecPrimZvPtEtaCent->GetAxis(2)->SetTitle("MC Eta");
 255   fMCRecPrimZvPtEtaCent->GetAxis(3)->SetTitle("Centrality");
 256   fMCRecPrimZvPtEtaCent->Sumw2();
 257
 258   fMCGenZvPtEtaCent = new THnSparseF("fMCGenZvPtEtaCent","mcZv:mcPt:mcEta:Centrality",4,binsZvPtEtaCent);
 259   fMCGenZvPtEtaCent->SetBinEdges(0,fBinsZv);
 260   fMCGenZvPtEtaCent->SetBinEdges(1,fBinsPt);
 261   fMCGenZvPtEtaCent->SetBinEdges(2,fBinsEta);
 262   fMCGenZvPtEtaCent->SetBinEdges(3,fBinsCentrality);
 263   fMCGenZvPtEtaCent->GetAxis(0)->SetTitle("MC Zv (cm)");
 264   fMCGenZvPtEtaCent->GetAxis(1)->SetTitle("MC Pt (GeV/c)");
 265   fMCGenZvPtEtaCent->GetAxis(2)->SetTitle("MC Eta");
 266   fMCGenZvPtEtaCent->GetAxis(3)->SetTitle("Centrality");
 267   fMCGenZvPtEtaCent->Sumw2();
 268
 269   fMCRecSecZvPtEtaCent = new THnSparseF("fMCRecSecZvPtEtaCent","mcZv:mcPt:mcEta:Centrality",4,binsZvPtEtaCent);
 270   fMCRecSecZvPtEtaCent->SetBinEdges(0,fBinsZv);
 271   fMCRecSecZvPtEtaCent->SetBinEdges(1,fBinsPt);
 272   fMCRecSecZvPtEtaCent->SetBinEdges(2,fBinsEta);
 273   fMCRecSecZvPtEtaCent->SetBinEdges(3,fBinsCentrality);
 274   fMCRecSecZvPtEtaCent->GetAxis(0)->SetTitle("MC Sec Zv (cm)");
 275   fMCRecSecZvPtEtaCent->GetAxis(1)->SetTitle("MC Sec Pt (GeV/c)");
 276   fMCRecSecZvPtEtaCent->GetAxis(2)->SetTitle("MC Sec Eta");
 277   fMCRecSecZvPtEtaCent->GetAxis(3)->SetTitle("Centrality");
 278   fMCRecSecZvPtEtaCent->Sumw2();
 279
 280   fMCRecPrimPhiPtEtaCent = new THnSparseF("fMCRecPrimPhiPtEtaCent","mcPhi:mcPt:mcEta:Centrality",4,binsPhiPtEtaCent);
 281   fMCRecPrimPhiPtEtaCent->SetBinEdges(0,fBinsPhi);
 282   fMCRecPrimPhiPtEtaCent->SetBinEdges(1,fBinsPt);
 283   fMCRecPrimPhiPtEtaCent->SetBinEdges(2,fBinsEta);
 284   fMCRecPrimPhiPtEtaCent->SetBinEdges(3,fBinsCentrality);
 285   fMCRecPrimPhiPtEtaCent->GetAxis(0)->SetTitle("MC Phi (cm)");
 286   fMCRecPrimPhiPtEtaCent->GetAxis(1)->SetTitle("MC Pt (GeV/c)");
 287   fMCRecPrimPhiPtEtaCent->GetAxis(2)->SetTitle("MC Eta");
 288   fMCRecPrimPhiPtEtaCent->GetAxis(3)->SetTitle("Centrality");
 289   fMCRecPrimPhiPtEtaCent->Sumw2();
 290
 291   fMCGenPhiPtEtaCent = new THnSparseF("fMCGenPhiPtEtaCent","mcPhi:mcPt:mcEta:Centrality",4,binsPhiPtEtaCent);
 292   fMCGenPhiPtEtaCent->SetBinEdges(0,fBinsPhi);
 293   fMCGenPhiPtEtaCent->SetBinEdges(1,fBinsPt);
 294   fMCGenPhiPtEtaCent->SetBinEdges(2,fBinsEta);
 295   fMCGenPhiPtEtaCent->SetBinEdges(3,fBinsCentrality);
 296   fMCGenPhiPtEtaCent->GetAxis(0)->SetTitle("MC Phi (cm)");
 297   fMCGenPhiPtEtaCent->GetAxis(1)->SetTitle("MC Pt (GeV/c)");
 298   fMCGenPhiPtEtaCent->GetAxis(2)->SetTitle("MC Eta");
 299   fMCGenPhiPtEtaCent->GetAxis(3)->SetTitle("Centrality");
 300   fMCGenPhiPtEtaCent->Sumw2();
 301
 302   fMCRecSecPhiPtEtaCent = new THnSparseF("fMCRecSecPhiPtEtaCent","mcPhi:mcPt:mcEta:Centrality",4,binsPhiPtEtaCent);
 303   fMCRecSecPhiPtEtaCent->SetBinEdges(0,fBinsPhi);
 304   fMCRecSecPhiPtEtaCent->SetBinEdges(1,fBinsPt);
 305   fMCRecSecPhiPtEtaCent->SetBinEdges(2,fBinsEta);
 306   fMCRecSecPhiPtEtaCent->SetBinEdges(3,fBinsCentrality);
 307   fMCRecSecPhiPtEtaCent->GetAxis(0)->SetTitle("MC Sec Phi (cm)");
 308   fMCRecSecPhiPtEtaCent->GetAxis(1)->SetTitle("MC Sec Pt (GeV/c)");
 309   fMCRecSecPhiPtEtaCent->GetAxis(2)->SetTitle("MC Sec Eta");
 310   fMCRecSecPhiPtEtaCent->GetAxis(3)->SetTitle("Centrality");
 311   fMCRecSecPhiPtEtaCent->Sumw2();
 312
 313   fPt = new TH1F("fPt","fPt",2000,0,200);
 314   fPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 315   fPt->GetYaxis()->SetTitle("dN/dp_{T}");
 316   fPt->Sumw2();
 317
 318   fMCPt = new TH1F("fMCPt","fMCPt",2000,0,200);
 319   fMCPt->GetXaxis()->SetTitle("MC p_{T} (GeV/c)");
 320   fMCPt->GetYaxis()->SetTitle("dN/dp_{T}");
 321   fMCPt->Sumw2();
 322
 323   fEventStatistics = new TH1F("fEventStatistics","fEventStatistics",10,0,10);
 324   fEventStatistics->GetYaxis()->SetTitle("number of events");
 325   fEventStatistics->SetBit(TH1::kCanRebin);
 326
 327   fEventStatisticsCentrality = new TH1F("fEventStatisticsCentrality","fEventStatisticsCentrality",fCentralityNbins-1, fBinsCentrality);
 328   fEventStatisticsCentrality->GetYaxis()->SetTitle("number of events");
 329
 330   fMCEventStatisticsCentrality = new TH1F("fMCEventStatisticsCentrality","fMCEventStatisticsCentrality",fCentralityNbins-1, fBinsCentrality);
 331   fMCEventStatisticsCentrality->GetYaxis()->SetTitle("number of MC events");
 332
 333   fAllEventStatisticsCentrality = new TH1F("fAllEventStatisticsCentrality","fAllEventStatisticsCentrality",fCentralityNbins-1, fBinsCentrality);
 334   fAllEventStatisticsCentrality->GetYaxis()->SetTitle("number of events");
 335
 336   fEventStatisticsCentralityTrigger = new TH2F("fEventStatisticsCentralityTrigger","fEventStatisticsCentralityTrigger;centrality;trigger",100,0,100,3,0,3);
 337   fEventStatisticsCentralityTrigger->Sumw2();
 338
 339   fZvMultCent = new THnSparseF("fZvMultCent","Zv:mult:Centrality",3,binsZvMultCent);
 340   fZvMultCent->SetBinEdges(0,fBinsZv);
 341   fZvMultCent->SetBinEdges(1,fBinsMult);
 342   fZvMultCent->SetBinEdges(2,fBinsCentrality);
 343   fZvMultCent->GetAxis(0)->SetTitle("Zv (cm)");
 344   fZvMultCent->GetAxis(1)->SetTitle("N_{acc}");
 345   fZvMultCent->GetAxis(2)->SetTitle("Centrality");
 346   fZvMultCent->Sumw2();
 347
 348   fTriggerStatistics = new TH1F("fTriggerStatistics","fTriggerStatistics",10,0,10);
 349   fTriggerStatistics->GetYaxis()->SetTitle("number of events");
 350
 351   fMCTrackPdgCode = new TH1F("fMCTrackPdgCode","fMCTrackPdgCode",100,0,10);
 352   fMCTrackPdgCode->GetYaxis()->SetTitle("number of tracks");
 353   fMCTrackPdgCode->SetBit(TH1::kCanRebin);
 354
 355   fMCTrackStatusCode = new TH1F("fMCTrackStatusCode","fMCTrackStatusCode",100,0,10);
 356   fMCTrackStatusCode->GetYaxis()->SetTitle("number of tracks");
 357   fMCTrackStatusCode->SetBit(TH1::kCanRebin);
 358
 359   fCharge = new TH1F("fCharge","fCharge",30, -5, 5);
 360   fCharge->GetXaxis()->SetTitle("Charge");
 361   fCharge->GetYaxis()->SetTitle("number of tracks");
 362
 363   fMCCharge = new TH1F("fMCCharge","fMCCharge",30, -5, 5);
 364   fMCCharge->GetXaxis()->SetTitle("MC Charge");
 365   fMCCharge->GetYaxis()->SetTitle("number of tracks");
 366
 367   fMCPdgPt = new TH2F("fMCPdgPt","fMCPdgPt",fPtNbins-1, fBinsPt, 100,0,100);
 368   fMCPdgPt->GetYaxis()->SetTitle("particle");
 369   fMCPdgPt->GetXaxis()->SetTitle("Pt (GeV/c)");
 370
 371   fMCHijingPrim = new TH1F("fMCHijingPrim","fMCHijingPrim",2,0,2);
 372   fMCHijingPrim->GetYaxis()->SetTitle("number of particles");
 373
 374
 375
 376   Int_t binsDCAxyDCAzPtEtaPhi[6] = { 200,200, fPtCheckNbins-1, fEtaCheckNbins-1, 18, fCentralityNbins-1};
 377   Double_t minDCAxyDCAzPtEtaPhi[6] = { -5, -5, 0, -1.5, 0., 0, };
 378   Double_t maxDCAxyDCAzPtEtaPhi[6] = { 5., 5., 100, 1.5, 2.*TMath::Pi(), 100};
 379
 380   fDCAPtAll = new THnSparseF("fDCAPtAll","fDCAPtAll",6, binsDCAxyDCAzPtEtaPhi, minDCAxyDCAzPtEtaPhi, maxDCAxyDCAzPtEtaPhi);
 381   fDCAPtAccepted = new THnSparseF("fDCAPtAccepted","fDCAPtAccepted",6, binsDCAxyDCAzPtEtaPhi, minDCAxyDCAzPtEtaPhi, maxDCAxyDCAzPtEtaPhi);
 382   fMCDCAPtSecondary = new THnSparseF("fMCDCAPtSecondary","fMCDCAPtSecondary",6, binsDCAxyDCAzPtEtaPhi, minDCAxyDCAzPtEtaPhi, maxDCAxyDCAzPtEtaPhi);
 383   fMCDCAPtPrimary = new THnSparseF("fMCDCAPtPrimary","fMCDCAPtPrimary",6, binsDCAxyDCAzPtEtaPhi, minDCAxyDCAzPtEtaPhi, maxDCAxyDCAzPtEtaPhi);
 384
 385   fDCAPtAll->SetBinEdges(2, fBinsPtCheck);
 386   fDCAPtAccepted->SetBinEdges(2, fBinsPtCheck);
 387   fMCDCAPtSecondary->SetBinEdges(2, fBinsPtCheck);
 388   fMCDCAPtPrimary->SetBinEdges(2, fBinsPtCheck);
 389
 390   fDCAPtAll->SetBinEdges(3, fBinsEtaCheck);
 391   fDCAPtAccepted->SetBinEdges(3, fBinsEtaCheck);
 392   fMCDCAPtSecondary->SetBinEdges(3, fBinsEtaCheck);
 393   fMCDCAPtPrimary->SetBinEdges(3, fBinsEtaCheck);
 394
 395   fDCAPtAll->SetBinEdges(5, fBinsCentrality);
 396   fDCAPtAccepted->SetBinEdges(5, fBinsCentrality);
 397   fMCDCAPtSecondary->SetBinEdges(5, fBinsCentrality);
 398   fMCDCAPtPrimary->SetBinEdges(5, fBinsCentrality);
 399
 400   fDCAPtAll->Sumw2();
 401   fDCAPtAccepted->Sumw2();
 402   fMCDCAPtSecondary->Sumw2();
 403   fMCDCAPtPrimary->Sumw2();
 404
 405   fDCAPtAll->GetAxis(0)->SetTitle("DCA_{xy} (cm)");
 406   fDCAPtAll->GetAxis(1)->SetTitle("DCA_{z} (cm)");
 407   fDCAPtAll->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
 408   fDCAPtAll->GetAxis(3)->SetTitle("#eta");
 409   fDCAPtAll->GetAxis(4)->SetTitle("#phi");
 410   fDCAPtAll->GetAxis(5)->SetTitle("Centrality");
 411
 412   fDCAPtAccepted->GetAxis(0)->SetTitle("DCA_{xy} (cm)");
 413   fDCAPtAccepted->GetAxis(1)->SetTitle("DCA_{z} (cm)");
 414   fDCAPtAccepted->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
 415   fDCAPtAccepted->GetAxis(3)->SetTitle("#eta");
 416   fDCAPtAccepted->GetAxis(4)->SetTitle("#phi");
 417   fDCAPtAccepted->GetAxis(5)->SetTitle("Centrality");
 418
 419   fMCDCAPtSecondary->GetAxis(0)->SetTitle("DCA_{xy} (cm)");
 420   fMCDCAPtSecondary->GetAxis(1)->SetTitle("DCA_{z} (cm)");
 421   fMCDCAPtSecondary->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
 422   fMCDCAPtSecondary->GetAxis(3)->SetTitle("#eta");
 423   fMCDCAPtSecondary->GetAxis(4)->SetTitle("#phi");
 424   fMCDCAPtSecondary->GetAxis(5)->SetTitle("Centrality");
 425
 426   fMCDCAPtPrimary->GetAxis(0)->SetTitle("DCA_{xy} (cm)");
 427   fMCDCAPtPrimary->GetAxis(1)->SetTitle("DCA_{z} (cm)");
 428   fMCDCAPtPrimary->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
 429   fMCDCAPtPrimary->GetAxis(3)->SetTitle("#eta");
 430   fMCDCAPtPrimary->GetAxis(4)->SetTitle("#phi");
 431   fMCDCAPtPrimary->GetAxis(5)->SetTitle("Centrality");
 432
 433
 434   char cFullTempTitle[255];
 435   char cTempTitleAxis0All[255];
 436   char cTempTitleAxis0Acc[255];
 437   //   char cTempTitleAxis1[255];
 438   char cFullTempName[255];
 439   char cTempNameAxis0[255];
 440   //   char cTempNameAxis1[255];
 441   const Int_t iNbinRowsClusters = 21;
 442   //   Double_t dBinsRowsClusters[iNbinRowsClusters] = {0, 7.95, 15.9, 23.85, 31.8, 39.75, 47.7, 55.65, 63.6, 71.55, 79.5, 87.45, 95.4, 103.35, 111.3, 119.25, 127.2, 135.15, 143.1, 151.05, 159.};
 443
 444   const Int_t iNbinChi = 51;
 445   const Int_t iNbinLength = 165;
 446   //   Double_t dBinsChi[iNbinChi] = {0, 0.2, 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2, 2.2, 2.4, 2.6, 2.8, 3, 3.2, 3.4, 3.6, 3.8, 4, 4.2, 4.4, 4.6, 4.8, 5, 5.2, 5.4, 5.6, 5.8, 6, 6.2, 6.4, 6.6, 6.8, 7, 7.2, 7.4, 7.6, 7.8, 8, 8.2, 8.4, 8.6, 8.8, 9, 9.2, 9.4, 9.6, 9.8,10.};
 447
 448   Int_t iNbin = 0;
 449   //   Double_t *dBins = 0x0;
 450   Double_t dBinMin = 0;
 451   Double_t dBinMax = 0;
 452
 453   for(Int_t iCheckQuant = 0; iCheckQuant < cqMax; iCheckQuant++)
 454   {
 455     // iCheckQuant: 0 = CrossedRows, 1 = Nclusters, 2 = Chi^2/clusterTPC
 456     if(iCheckQuant == cqCrossedRows)
 457     {
 458       snprintf(cTempTitleAxis0All,255, "NcrossedRows before Cut");
 459       snprintf(cTempTitleAxis0Acc,255, "NcrossedRows after Cut");
 460       snprintf(cTempNameAxis0,255, "CrossedRows");
 461       iNbin = iNbinRowsClusters;
 462       dBinMin = 0;
 463       dBinMax = 159.;
 464     }
 465     else if(iCheckQuant == cqNcluster)
 466     {
 467       snprintf(cTempTitleAxis0All,255, "Nclusters before Cut");
 468       snprintf(cTempTitleAxis0Acc,255, "Nclusters after Cut");
 469       snprintf(cTempNameAxis0,255, "Clusters");
 470       iNbin = iNbinRowsClusters;
 471       dBinMin = 0;
 472       dBinMax = 159.;
 473     }
 474     else if(iCheckQuant == cqChi)
 475     {
 476       snprintf(cTempTitleAxis0All,255, "#Chi^{2}/cluster before Cut");
 477       snprintf(cTempTitleAxis0Acc,255, "#Chi^{2}/cluster after Cut");
 478       snprintf(cTempNameAxis0,255, "Chi");
 479       iNbin = iNbinChi;
 480       dBinMin = 0;
 481       dBinMax = 10.;
 482     }
 483     else if(iCheckQuant == cqLength)
 484     {
 485       snprintf(cTempTitleAxis0All,255, "Length in TPC before Cut (cm)");
 486       snprintf(cTempTitleAxis0Acc,255, "Length in TPC after Cut (cm)");
 487       snprintf(cTempNameAxis0,255, "Length");
 488       iNbin = iNbinLength;
 489       dBinMin = 0;
 490       dBinMax = 165.;
 491     }
 492
 493     Int_t binsCheckPtEtaPhi[5] = { iNbin, fPtCheckNbins-1, fEtaCheckNbins-1, 18, fCentralityNbins-1};
 494 //     Int_t binsCheckPtEtaPhi[5] = { iNbin, fPtNbins-1, fEtaCheckNbins-1, 18, fCentralityNbins-1};
 495     Double_t minCheckPtEtaPhi[5] = { dBinMin,  0, -1.5, 0., 0, };
 496     Double_t maxCheckPtEtaPhi[5] = { dBinMax, 100, 1.5, 2.*TMath::Pi(), 100};
 497
 498     snprintf(cFullTempName, 255, "f%sPtEtaPhiAll",cTempNameAxis0);
 499     snprintf(cFullTempTitle, 255,"%s;%s;p_{T} (GeV/c);#eta;#phi;Centrality", cFullTempName, cTempTitleAxis0All);
 500     fCrossCheckAll[iCheckQuant] = new THnF(cFullTempName, cFullTempTitle, 5, binsCheckPtEtaPhi, minCheckPtEtaPhi, maxCheckPtEtaPhi);
 501     fCrossCheckAll[iCheckQuant]->SetBinEdges(1, fBinsPtCheck);
 502     fCrossCheckAll[iCheckQuant]->SetBinEdges(2, fBinsEtaCheck);
 503     fCrossCheckAll[iCheckQuant]->Sumw2();
 504
 505     snprintf(cFullTempName, 255, "f%sPtEtaPhiAcc",cTempNameAxis0);
 506     snprintf(cFullTempTitle, 255,"%s;%s;p_{T} (GeV/c);#eta;#phi;Centrality", cFullTempName, cTempTitleAxis0Acc);
 507     fCrossCheckAcc[iCheckQuant] = new THnF(cFullTempName, cFullTempTitle, 5, binsCheckPtEtaPhi, minCheckPtEtaPhi, maxCheckPtEtaPhi);
 508     fCrossCheckAcc[iCheckQuant]->SetBinEdges(1, fBinsPtCheck);
 509     fCrossCheckAcc[iCheckQuant]->SetBinEdges(2, fBinsEtaCheck);
 510     fCrossCheckAcc[iCheckQuant]->Sumw2();
 511   } // end iCheckQuant
 512
 513   fCutPercClusters = new TH1F("fCutPercClusters","fCutPercClusters;NclustersTPC;counts",160,0,160);
 514   fCutPercClusters->Sumw2();
 515   fCutPercCrossed = new TH1F("fCutPercCrossed","fCutPercCrossed;NcrossedRowsTPC;counts",160,0,160);
 516   fCutPercCrossed->Sumw2();
 517
 518   fCrossCheckRowsLength = new TH2F("fCrossCheckRowsLength","fCrossCheckRowsLength;Length in TPC;NcrossedRows",170,0,170,170,0,170);
 519   fCrossCheckRowsLength->Sumw2();
 520
 521   fCrossCheckClusterLength = new TH2F("fCrossCheckClusterLength","fCrossCheckClusterLength;Length in TPC;NClusters",170,0,170,170,0,170);
 522   fCrossCheckClusterLength->Sumw2();
 523
 524   fCrossCheckRowsLengthAcc = new TH2F("fCrossCheckRowsLengthAcc","fCrossCheckRowsLengthAcc;Length in TPC;NcrossedRows",170,0,170,170,0,170);
 525   fCrossCheckRowsLengthAcc->Sumw2();
 526
 527   fCrossCheckClusterLengthAcc = new TH2F("fCrossCheckClusterLengthAcc","fCrossCheckClusterLengthAcc;Length in TPC;NClusters",170,0,170,170,0,170);
 528   fCrossCheckClusterLengthAcc->Sumw2();
 529
 530
 531   // Add Histos, Profiles etc to List
 532   fOutputList->Add(fZvPtEtaCent);
 533   fOutputList->Add(fPhiPtEtaCent);
 534   fOutputList->Add(fPt);
 535   fOutputList->Add(fMCRecPrimZvPtEtaCent);
 536   fOutputList->Add(fMCGenZvPtEtaCent);
 537   fOutputList->Add(fMCRecSecZvPtEtaCent);
 538   fOutputList->Add(fMCRecPrimPhiPtEtaCent);
 539   fOutputList->Add(fMCGenPhiPtEtaCent);
 540   fOutputList->Add(fMCRecSecPhiPtEtaCent);
 541   fOutputList->Add(fMCPt);
 542   fOutputList->Add(fEventStatistics);
 543   fOutputList->Add(fEventStatisticsCentrality);
 544   fOutputList->Add(fMCEventStatisticsCentrality);
 545   fOutputList->Add(fAllEventStatisticsCentrality);
 546   fOutputList->Add(fEventStatisticsCentralityTrigger);
 547   fOutputList->Add(fZvMultCent);
 548   fOutputList->Add(fTriggerStatistics);
 549   fOutputList->Add(fMCTrackPdgCode);
 550   fOutputList->Add(fMCTrackStatusCode);
 551   fOutputList->Add(fCharge);
 552   fOutputList->Add(fMCCharge);
 553   fOutputList->Add(fMCPdgPt);
 554   fOutputList->Add(fMCHijingPrim);
 555   fOutputList->Add(fDCAPtAll);
 556   fOutputList->Add(fDCAPtAccepted);
 557   fOutputList->Add(fMCDCAPtSecondary);
 558   fOutputList->Add(fMCDCAPtPrimary);
 559   for(Int_t i = 0; i < cqMax; i++)
 560   {
 561     fOutputList->Add(fCrossCheckAll[i]);
 562     fOutputList->Add(fCrossCheckAcc[i]);
 563   }
 564   fOutputList->Add(fCutPercClusters);
 565   fOutputList->Add(fCutPercCrossed);
 566   fOutputList->Add(fCrossCheckRowsLength);
 567   fOutputList->Add(fCrossCheckClusterLength);
 568   fOutputList->Add(fCrossCheckRowsLengthAcc);
 569   fOutputList->Add(fCrossCheckClusterLengthAcc);
 570
 571   PostData(1, fOutputList);
 572 }
 573
 574 void AlidNdPtAnalysisPbPbAOD::UserExec(Option_t *option)
 575 {
 576
 577   // Main Loop
 578   // called for each event
 579   fEventStatistics->Fill("all events",1);
 580
 581   // set ZERO pointers:
 582   AliInputEventHandler *inputHandler = NULL;
 583   AliAODTrack *track = NULL;
 584   AliAODMCParticle *mcPart = NULL;
 585   AliAODMCHeader *mcHdr = NULL;
 586   AliGenHijingEventHeader *genHijingHeader = NULL;
 587   //AliGenPythiaEventHeader *genPythiaHeader = NULL;
 588
 589   Bool_t bIsEventSelectedMB = kFALSE;
 590   Bool_t bIsEventSelectedSemi = kFALSE;
 591   Bool_t bIsEventSelectedCentral = kFALSE;
 592   Bool_t bIsEventSelected = kFALSE;
 593   Bool_t bIsPrimary = kFALSE;
 594   Bool_t bIsHijingParticle = kFALSE;
 595   Bool_t bMotherIsHijingParticle = kFALSE;
 596   //Bool_t bIsPythiaParticle = kFALSE;
 597   Bool_t bEventHasATrack = kFALSE;
 598   Bool_t bEventHasATrackInRange = kFALSE;
 599   Int_t nTriggerFired = 0;
 600
 601
 602   Double_t dMCTrackZvPtEtaCent[4] = {0};
 603   Double_t dTrackZvPtEtaCent[4] = {0};
 604
 605   Double_t dMCTrackPhiPtEtaCent[4] = {0};
 606   Double_t dTrackPhiPtEtaCent[4] = {0};
 607
 608   Double_t dDCA[2] = {0};
 609
 610   Double_t dMCEventZv = -100;
 611   Double_t dEventZv = -100;
 612   Int_t iAcceptedMultiplicity = 0;
 613
 614   fIsMonteCarlo = kFALSE;
 615
 616   AliAODEvent *eventAOD = 0x0;
 617   eventAOD = dynamic_cast<AliAODEvent*>( InputEvent() );
 618   if (!eventAOD) {
 619     AliWarning("ERROR: eventAOD not available \n");
 620     return;
 621   }
 622
 623   // check, which trigger has been fired
 624   inputHandler = (AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 625   bIsEventSelectedMB = ( inputHandler->IsEventSelected() & AliVEvent::kMB);
 626   bIsEventSelectedSemi = ( inputHandler->IsEventSelected() & AliVEvent::kSemiCentral);
 627   bIsEventSelectedCentral = ( inputHandler->IsEventSelected() & AliVEvent::kCentral);
 628
 629   if(bIsEventSelectedMB || bIsEventSelectedSemi || bIsEventSelectedCentral) fTriggerStatistics->Fill("all triggered events",1);
 630   if(bIsEventSelectedMB) { fTriggerStatistics->Fill("MB trigger",1); nTriggerFired++; }
 631   if(bIsEventSelectedSemi) { fTriggerStatistics->Fill("SemiCentral trigger",1); nTriggerFired++; }
 632   if(bIsEventSelectedCentral) { fTriggerStatistics->Fill("Central trigger",1); nTriggerFired++; }
 633   if(nTriggerFired == 0) { fTriggerStatistics->Fill("No trigger",1); }
 634
 635   bIsEventSelected = ( inputHandler->IsEventSelected() & GetCollisionCandidates() );
 636
 637   // only take tracks of events, which are triggered
 638   if(nTriggerFired == 0) { return; }
 639
 640   //   if( !bIsEventSelected || nTriggerFired>1 ) return;
 641
 642   //   fEventStatistics->Fill("events with only coll. cand.", 1);
 643
 644
 645
 646   // check if there is a stack, if yes, then do MC loop
 647   TList *list = eventAOD->GetList();
 648   TClonesArray *stack = 0x0;
 649   stack = (TClonesArray*)list->FindObject(AliAODMCParticle::StdBranchName());
 650
 651   if( stack )
 652   {
 653     fIsMonteCarlo = kTRUE;
 654
 655     mcHdr = (AliAODMCHeader*)list->FindObject(AliAODMCHeader::StdBranchName());
 656
 657     genHijingHeader = GetHijingEventHeader(mcHdr);
 658     //     genPythiaHeader = GetPythiaEventHeader(mcHdr);
 659
 660     if(!genHijingHeader) { return; }
 661
 662     //     if(!genPythiaHeader)  { return; }
 663
 664     dMCEventZv = mcHdr->GetVtxZ();
 665     dMCTrackZvPtEtaCent[0] = dMCEventZv;
 666     fEventStatistics->Fill("MC all events",1);
 667   }
 668
 669   AliCentrality* aCentrality = eventAOD->GetCentrality();
 670   Double_t dCentrality = aCentrality->GetCentralityPercentile("V0M");
 671
 672   if( dCentrality < 0 ) return;
 673   fEventStatistics->Fill("after centrality selection",1);
 674
 675
 676
 677   // start with MC truth analysis
 678   if(fIsMonteCarlo)
 679   {
 680
 681     if( dMCEventZv > GetCutMaxZVertex() )  { return; }
 682
 683     dMCTrackZvPtEtaCent[0] = dMCEventZv;
 684
 685     fEventStatistics->Fill("MC afterZv cut",1);
 686
 687     for(Int_t iMCtrack = 0; iMCtrack < stack->GetEntriesFast(); iMCtrack++)
 688     {
 689       mcPart =(AliAODMCParticle*)stack->At(iMCtrack);
 690
 691       // check for charge
 692       if( !(IsMCTrackAccepted(mcPart)) ) continue;
 693
 694       if(!IsHijingParticle(mcPart, genHijingHeader)) { continue; }
 695
 696       if(mcPart->IsPhysicalPrimary() )
 697       {
 698         fMCHijingPrim->Fill("IsPhysicalPrimary",1);
 699       }
 700       else
 701       {
 702         fMCHijingPrim->Fill("NOT a primary",1);
 703         continue;
 704       }
 705
 706
 707       //
 708       // ======================== fill histograms ========================
 709       dMCTrackZvPtEtaCent[1] = mcPart->Pt();
 710       dMCTrackZvPtEtaCent[2] = mcPart->Eta();
 711       dMCTrackZvPtEtaCent[3] = dCentrality;
 712       fMCGenZvPtEtaCent->Fill(dMCTrackZvPtEtaCent);
 713
 714       dMCTrackPhiPtEtaCent[0] = mcPart->Phi();
 715       dMCTrackPhiPtEtaCent[1] = mcPart->Pt();
 716       dMCTrackPhiPtEtaCent[2] = mcPart->Eta();
 717       dMCTrackPhiPtEtaCent[3] = dCentrality;
 718       fMCGenPhiPtEtaCent->Fill(dMCTrackPhiPtEtaCent);
 719
 720       bEventHasATrack = kTRUE;
 721
 722
 723       if( (dMCTrackZvPtEtaCent[1] > GetCutPtMin() ) &&
 724         (dMCTrackZvPtEtaCent[1] < GetCutPtMax() ) &&
 725         (dMCTrackZvPtEtaCent[2] > GetCutEtaMin() ) &&
 726         (dMCTrackZvPtEtaCent[2] < GetCutEtaMax() ) )
 727       {
 728         fMCPt->Fill(mcPart->Pt());
 729         fMCCharge->Fill(mcPart->Charge()/3.);
 730         bEventHasATrackInRange = kTRUE;
 731       }
 732
 733     }
 734   } // isMonteCarlo
 735
 736   if(bEventHasATrack) { fEventStatistics->Fill("MC events with tracks",1); }
 737   if(bEventHasATrackInRange)
 738   {
 739     fEventStatistics->Fill("MC events with tracks in range",1);
 740     fMCEventStatisticsCentrality->Fill(dCentrality);
 741   }
 742   bEventHasATrack = kFALSE;
 743   bEventHasATrackInRange = kFALSE;
 744
 745
 746
 747   // Loop over recontructed tracks
 748
 749   dEventZv = eventAOD->GetPrimaryVertex()->GetZ();
 750   if( TMath::Abs(dEventZv) > GetCutMaxZVertex() ) return;
 751
 752   fAllEventStatisticsCentrality->Fill(dCentrality/*, nTriggerFired*/);
 753
 754   fEventStatistics->Fill("after Zv cut",1);
 755
 756   dTrackZvPtEtaCent[0] = dEventZv;
 757
 758   if(AreRelativeCutsEnabled())
 759   {
 760     if(!SetRelativeCuts(eventAOD)) return;
 761   }
 762
 763   for(Int_t itrack = 0; itrack < eventAOD->GetNumberOfTracks(); itrack++)
 764   {
 765     track = eventAOD->GetTrack(itrack);
 766     if(!track) continue;
 767
 768     mcPart = NULL;
 769     dMCTrackZvPtEtaCent[1] = 0;
 770     dMCTrackZvPtEtaCent[2] = 0;
 771     dMCTrackZvPtEtaCent[3] = 0;
 772
 773     dMCTrackPhiPtEtaCent[0] = 0;
 774     dMCTrackPhiPtEtaCent[1] = 0;
 775     dMCTrackPhiPtEtaCent[2] = 0;
 776     dMCTrackPhiPtEtaCent[3] = 0;
 777
 778     bIsPrimary = kFALSE;
 779
 780     GetDCA(track, eventAOD, dDCA);
 781
 782     Double_t dDCAxyDCAzPt[5] = { dDCA[0], dDCA[1], track->Pt(), track->Eta(), track->Phi() };
 783
 784     fDCAPtAll->Fill(dDCAxyDCAzPt);
 785
 786     if( !(IsTrackAccepted(track, dCentrality, eventAOD->GetMagneticField())) ) continue;
 787
 788     dTrackZvPtEtaCent[1] = track->Pt();
 789     dTrackZvPtEtaCent[2] = track->Eta();
 790     dTrackZvPtEtaCent[3] = dCentrality;
 791
 792     dTrackPhiPtEtaCent[0] = track->Phi();
 793     dTrackPhiPtEtaCent[1] = track->Pt();
 794     dTrackPhiPtEtaCent[2] = track->Eta();
 795     dTrackPhiPtEtaCent[3] = dCentrality;
 796
 797     if( fIsMonteCarlo )
 798     {
 799       mcPart = (AliAODMCParticle*)stack->At(TMath::Abs(track->GetLabel()));
 800       if( !mcPart ) { continue; }
 801
 802       // check for charge
 803       // if( !(IsMCTrackAccepted(mcPart)) ) {  continue; }
 804
 805       bIsHijingParticle = IsHijingParticle(mcPart, genHijingHeader);
 806       //       bIsPythiaParticle = IsPythiaParticle(mcPart, genPythiaHeader);
 807
 808       //       if(!bIsHijingParticle) continue; // only take real tracks, not injected ones
 809
 810       bIsPrimary = mcPart->IsPhysicalPrimary();
 811
 812       dMCTrackZvPtEtaCent[1] = mcPart->Pt();
 813       dMCTrackZvPtEtaCent[2] = mcPart->Eta();
 814       dMCTrackZvPtEtaCent[3] = dCentrality;
 815
 816       dMCTrackPhiPtEtaCent[0] = mcPart->Phi();
 817     dMCTrackPhiPtEtaCent[1] = mcPart->Pt();
 818     dMCTrackPhiPtEtaCent[2] = mcPart->Eta();
 819     dMCTrackPhiPtEtaCent[3] = dCentrality;
 820
 821       if(bIsPrimary && bIsHijingParticle)
 822       {
 823         fMCRecPrimZvPtEtaCent->Fill(dMCTrackZvPtEtaCent);
 824         fMCRecPrimPhiPtEtaCent->Fill(dMCTrackPhiPtEtaCent);
 825         if( (TMath::Abs(mcPart->Eta()) < 0.8) && (dCentrality<5.) ) { fMCPdgPt->Fill(mcPart->Pt(), Form("%s",GetParticleName(mcPart->GetPdgCode())), 1); }
 826         fMCDCAPtPrimary->Fill(dDCAxyDCAzPt);
 827       }
 828
 829       if(!bIsPrimary /*&& !bIsHijingParticle*/)
 830       {
 831         Int_t indexMoth = mcPart->GetMother();
 832         if(indexMoth >= 0)
 833         {
 834           AliAODMCParticle* moth = (AliAODMCParticle*)stack->At(indexMoth);
 835           bMotherIsHijingParticle = IsHijingParticle(moth, genHijingHeader);
 836
 837           if(bMotherIsHijingParticle) // only store secondaries, which come from a not embedded signal!
 838           {
 839             fMCTrackStatusCode->Fill(Form("%d",mcPart->GetStatus()), 1);
 840
 841
 842             fMCRecSecZvPtEtaCent->Fill(dMCTrackZvPtEtaCent);
 843             fMCRecSecPhiPtEtaCent->Fill(dMCTrackPhiPtEtaCent);
 844             fMCDCAPtSecondary->Fill(dDCAxyDCAzPt);
 845             fMCTrackPdgCode->Fill(Form("%s_H%i_H%i",GetParticleName(moth->GetPdgCode()),bMotherIsHijingParticle, bIsHijingParticle), 1);
 846             //    delete moth;
 847           }
 848         }
 849       }
 850     } // end isMonteCarlo
 851
 852     // ======================== fill histograms ========================
 853
 854     // only keep prim and sec from not embedded signal
 855     Bool_t bKeepMCTrack = kFALSE;
 856     if(fIsMonteCarlo)
 857     {
 858       if( (bIsHijingParticle && bIsPrimary) ^ (bMotherIsHijingParticle && !bIsPrimary) )
 859       {
 860         bKeepMCTrack = kTRUE;
 861       }
 862       else
 863       {
 864         continue;
 865       }
 866     }
 867
 868     bEventHasATrack = kTRUE;
 869
 870     fZvPtEtaCent->Fill(dTrackZvPtEtaCent);
 871     fPhiPtEtaCent->Fill(dTrackPhiPtEtaCent);
 872
 873     fDCAPtAccepted->Fill(dDCAxyDCAzPt);
 874
 875     if( (dTrackZvPtEtaCent[1] > GetCutPtMin()) &&
 876       (dTrackZvPtEtaCent[1] < GetCutPtMax()) &&
 877       (dTrackZvPtEtaCent[2] > GetCutEtaMin()) &&
 878       (dTrackZvPtEtaCent[2] < GetCutEtaMax()) )
 879     {
 880       iAcceptedMultiplicity++;
 881       bEventHasATrackInRange = kTRUE;
 882       fPt->Fill(track->Pt());
 883       fCharge->Fill(track->Charge());
 884     }
 885   } // end track loop
 886
 887   if(bEventHasATrack) { fEventStatistics->Fill("events with tracks",1); bEventHasATrack = kFALSE; }
 888
 889   if(bEventHasATrackInRange)
 890   {
 891     fEventStatistics->Fill("events with tracks in range",1);
 892     fEventStatisticsCentrality->Fill(dCentrality);
 893     bEventHasATrackInRange = kFALSE;
 894
 895     if(bIsEventSelectedMB) fEventStatisticsCentralityTrigger->Fill(dCentrality, 0);
 896     if(bIsEventSelectedSemi) fEventStatisticsCentralityTrigger->Fill(dCentrality, 1);
 897     if(bIsEventSelectedCentral) fEventStatisticsCentralityTrigger->Fill(dCentrality, 2);
 898   }
 899
 900   Double_t dEventZvMultCent[3] = {dEventZv, iAcceptedMultiplicity, dCentrality};
 901   fZvMultCent->Fill(dEventZvMultCent);
 902
 903   PostData(1, fOutputList);
 904
 905 }
 906
 907 Bool_t AlidNdPtAnalysisPbPbAOD::SetRelativeCuts(AliAODEvent *event)
 908 {
 909   if(!event) return kFALSE;
 910
 911   AliAODTrack *tr = 0x0;
 912   TH1F *hCluster = new TH1F("hCluster","hCluster",160,0,160);
 913   TH1F *hCrossed = new TH1F("hCrossed","hCrossed",160,0,160);
 914
 915   for(Int_t itrack = 0; itrack < event->GetNumberOfTracks(); itrack++)
 916   {
 917     tr = event->GetTrack(itrack);
 918     if(!tr) continue;
 919
 920     // do some selection already
 921     //if(!(tr->TestFilterBit(AliAODTrack::kTrkGlobal)) ) { continue; }
 922
 923     Double_t dNClustersTPC = tr->GetTPCNcls();
 924     Double_t dCrossedRowsTPC = tr->GetTPCClusterInfo(2,1);
 925
 926     hCluster->Fill(dNClustersTPC);
 927     hCrossed->Fill(dCrossedRowsTPC);
 928   }
 929
 930   // loop trough histogram to check, where percentage is reach
 931   Double_t dTotIntCluster = hCluster->Integral();
 932   Double_t dTotIntCrossed = hCrossed->Integral();
 933   Float_t dIntCluster = 0;
 934   Float_t dIntCrossed = 0;
 935
 936   if(dTotIntCluster)
 937   {
 938     for(Int_t i = 0; i < hCluster->GetNbinsX(); i++)
 939     {
 940       if(hCluster->GetBinCenter(i) < 0) continue;
 941       dIntCluster += hCluster->GetBinContent(i);
 942       if(dIntCluster/dTotIntCluster > (1-GetCutPercMinNClustersTPC()))
 943       {
 944         SetCutMinNClustersTPC(hCluster->GetBinCenter(i));
 945         fCutPercClusters->Fill(hCluster->GetBinCenter(i));
 946         break;
 947       }
 948     }
 949   }
 950
 951   if(dTotIntCrossed)
 952   {
 953     for(Int_t i = 0; i < hCrossed->GetNbinsX(); i++)
 954     {
 955       if(hCrossed->GetBinCenter(i) < 0) continue;
 956       dIntCrossed += hCrossed->GetBinContent(i);
 957       if(dIntCrossed/dTotIntCrossed > (1-GetCutPercMinNCrossedRowsTPC()))
 958       {
 959         SetCutMinNClustersTPC(hCrossed->GetBinCenter(i));
 960         fCutPercCrossed->Fill(hCrossed->GetBinCenter(i));
 961         break;
 962       }
 963     }
 964   }
 965
 966   delete hCrossed;
 967   delete hCluster;
 968   return kTRUE;
 969
 970 }
 971
 972 Bool_t AlidNdPtAnalysisPbPbAOD::IsTrackAccepted(AliAODTrack *tr, Double_t dCentrality, Double_t bMagZ)
 973 {
 974   if(!tr) return kFALSE;
 975
 976   if(tr->Charge()==0) { return kFALSE; }
 977
 978   //
 979   // as done in AliAnalysisTaskFragmentationFunction
 980   //
 981
 982   Short_t sign = tr->Charge();
 983   Double_t xyz[50];
 984   Double_t pxpypz[50];
 985   Double_t cv[100];
 986
 987   for(Int_t i = 0; i < 100; i++) cv[i] = 0;
 988   for(Int_t i = 0; i < 50; i++) xyz[i] = 0;
 989   for(Int_t i = 0; i < 50; i++) pxpypz[i] = 0;
 990
 991   tr->GetXYZ(xyz);
 992   tr->GetPxPyPz(pxpypz);
 993
 994   // similar error occured as this one:
 995   // See https://savannah.cern.ch/bugs/?102721
 996   // which is one of the two 11h re-filtering follow-ups:
 997   // Andrea Dainese now first does the beam pipe
 998   // check and then copies from the vtrack (was the other
 999   // way around) to avoid the crash in the etp::Set()
1000
1001 //   if(xyz[0]*xyz[0]+xyz[1]*xyz[1] > 3.*3.) { return kFALSE; }
1002
1003   AliExternalTrackParam * par = new AliExternalTrackParam(xyz, pxpypz, cv, sign);
1004   if(!par) { return kFALSE; }
1005   AliESDtrack dummy;
1006 //   Double_t dLength = dummy.GetLengthInActiveZone(par,3,236, -5 ,0,0);
1007 //   Double_t dLengthInTPC = GetLengthInTPC(tr, 1.8, 220, bMagZ);
1008
1009   Double_t dLengthInTPC = dummy.GetLengthInActiveZone(par,3,236, bMagZ ,0,0);
1010   Double_t dNClustersTPC = tr->GetTPCNcls();
1011   Double_t dCrossedRowsTPC = tr->GetTPCClusterInfo(2,1);
1012   Double_t dFindableClustersTPC = tr->GetTPCNclsF();
1013   Double_t dChi2PerClusterTPC = (dNClustersTPC>0)?tr->Chi2perNDF()*(dNClustersTPC-5)/dNClustersTPC:-1.; // see AliDielectronVarManager.h
1014
1015   //   hAllCrossedRowsTPC->Fill(dCrossedRowsTPC);
1016
1017   Double_t dCheck[cqMax] = {dCrossedRowsTPC, dNClustersTPC, dChi2PerClusterTPC, dLengthInTPC};// = new Double_t[cqMax];
1018   Double_t dKine[kqMax] = {tr->Pt(), tr->Eta(), tr->Phi()};// = new Double_t[kqMax];
1019   //   dKine[0] = tr->Pt();
1020   //   dKine[1] = tr->Eta();
1021   //   dKine[2] = tr->Phi();
1022   //
1023   //   dCheck[0] = dCrossedRowsTPC;
1024   //   dCheck[1] = dNClustersTPC;
1025   //   dCheck[2] = dChi2PerClusterTPC;
1026
1027
1028   FillDebugHisto(dCheck, dKine, dCentrality, kFALSE);
1029
1030   // first cut on length
1031
1032   if( DoCutLengthInTPCPtDependent() && ( dLengthInTPC < GetPrefactorLengthInTPCPtDependent()*(130-5*TMath::Abs(1./tr->Pt())) )  ) { return kFALSE; }
1033
1034   // filter bit 5
1035   if(!(tr->TestFilterBit(AliAODTrack::kTrkGlobal)) ) { return kFALSE; }
1036
1037   // filter bit 4
1038   //   if(!(tr->TestFilterBit(AliAODTrack::kTrkGlobalNoDCA)) ) { return kFALSE; }
1039
1040   //   hFilterCrossedRowsTPC->Fill(dCrossedRowsTPC);
1041
1042
1043    if(dFindableClustersTPC == 0) {return kFALSE; }
1044    if(dCrossedRowsTPC < GetCutMinNCrossedRowsTPC()) { return kFALSE; }
1045    if( (dCrossedRowsTPC/dFindableClustersTPC) < GetCutMinRatioCrossedRowsOverFindableClustersTPC() ) { return kFALSE; }
1046    if(dNClustersTPC < GetCutMinNClustersTPC()) { return kFALSE; }
1047
1048   if (!(tr->GetStatus() & AliVTrack::kITSrefit)) { return kFALSE; } // no ITS refit
1049
1050   // do a relativ cut in Nclusters, both time at 80% of mean
1051   //   if(fIsMonteCarlo)
1052   //   {
1053     //     if(dNClustersTPC < 88) { return kFALSE; }
1054     //   }
1055     //   else
1056     //   {
1057       //     if(dNClustersTPC < 76) { return kFALSE; }
1058       //   }
1059
1060
1061       FillDebugHisto(dCheck, dKine, dCentrality, kTRUE);
1062
1063
1064       //   hAccNclsTPC->Fill(dNClustersTPC);
1065       //   hAccCrossedRowsTPC->Fill(dCrossedRowsTPC);
1066       //   Double_t dFindableClustersTPC = tr->GetTPCNclsF();
1067       //   Double_t dChi2PerClusterTPC = (dNClustersTPC>0)?tr->Chi2perNDF()*(dNClustersTPC-5)/dNClustersTPC:-1.; // see AliDielectronVarManager.h
1068       //
1069       //   Bool_t bIsFromKink = kFALSE;
1070       //   if(tr->GetProdVertex()->GetType() == AliAODVertex::kKink) bIsFromKink = kTRUE;
1071       //
1072       //   // from AliAnalysisTaskPIDqa.cxx
1073       //   ULong_t uStatus = tr->GetStatus();
1074       //   Bool_t bHasRefitTPC = kFALSE;
1075       //   Bool_t bHasRefitITS = kFALSE;
1076       //
1077       //   if ((uStatus & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) bHasRefitTPC = kTRUE;
1078       //   if ((uStatus & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) bHasRefitITS = kTRUE;
1079       //
1080       //   // from AliDielectronVarManager.h
1081       //   Bool_t bHasHitInSPD = kFALSE;
1082       //   for (Int_t iC=0; iC<2; iC++)
1083       //   {
1084         //     if (((tr->GetITSClusterMap()) & (1<<(iC))) > 0) {  bHasHitInSPD = kTRUE;  }
1085         //   }
1086         //
1087         //   Double_t dNClustersITS = tr->GetITSNcls();
1088
1089         // cuts to be done:
1090         // TPC
1091         //   esdTrackCuts->SetMinNCrossedRowsTPC(70);
1092         //   esdTrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
1093         //
1094         //   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
1095         //   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
1096         //   esdTrackCuts->SetRequireTPCRefit(kTRUE);
1097         // ITS
1098         //   esdTrackCuts->SetRequireITSRefit(kTRUE);
1099         //   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,       AliESDtrackCuts::kAny);
1100         //
1101         //   esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0105+0.0350/pt^1.1");
1102         //   esdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
1103         //
1104         //   esdTrackCuts->SetMaxDCAToVertexZ(2);
1105         //   esdTrackCuts->SetDCAToVertex2D(kFALSE);
1106         //   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1107         //
1108         //   esdTrackCuts->SetMaxChi2PerClusterITS(36);
1109
1110         //     delete [] dKine;
1111         //     delete [] dCheck;
1112         return kTRUE;
1113 }
1114
1115 Bool_t AlidNdPtAnalysisPbPbAOD::FillDebugHisto(Double_t *dCrossCheckVar, Double_t *dKineVar, Double_t dCentrality, Bool_t bIsAccepted)
1116 {
1117   if(bIsAccepted)
1118   {
1119     for(Int_t iCrossCheck = 0; iCrossCheck < cqMax; iCrossCheck++)
1120     {
1121       Double_t dFillIt[5] = {dCrossCheckVar[iCrossCheck], dKineVar[0], dKineVar[1], dKineVar[2], dCentrality};
1122       fCrossCheckAcc[iCrossCheck]->Fill(dFillIt);
1123     }
1124
1125     fCrossCheckRowsLengthAcc->Fill(dCrossCheckVar[cqLength], dCrossCheckVar[cqCrossedRows]);
1126     fCrossCheckClusterLengthAcc->Fill(dCrossCheckVar[cqLength], dCrossCheckVar[cqNcluster]);
1127   }
1128   else
1129   {
1130     for(Int_t iCrossCheck = 0; iCrossCheck < cqMax; iCrossCheck++)
1131     {
1132       Double_t dFillIt[5] = {dCrossCheckVar[iCrossCheck], dKineVar[0], dKineVar[1], dKineVar[2], dCentrality};
1133       fCrossCheckAll[iCrossCheck]->Fill(dFillIt);
1134     }
1135
1136     fCrossCheckRowsLength->Fill(dCrossCheckVar[cqLength], dCrossCheckVar[cqCrossedRows]);
1137     fCrossCheckClusterLength->Fill(dCrossCheckVar[cqLength], dCrossCheckVar[cqNcluster]);
1138   }
1139
1140   return kTRUE;
1141
1142 }
1143
1144 Bool_t AlidNdPtAnalysisPbPbAOD::GetDCA(const AliAODTrack *track, AliAODEvent *evt, Double_t d0z0[2])
1145 {
1146   // function adapted from AliDielectronVarManager.h
1147
1148   if(track->TestBit(AliAODTrack::kIsDCA)){
1149     d0z0[0]=track->DCA();
1150     d0z0[1]=track->ZAtDCA();
1151     return kTRUE;
1152   }
1153
1154   Bool_t ok=kFALSE;
1155   if(evt) {
1156     Double_t covd0z0[3];
1157     //AliAODTrack copy(*track);
1158     AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1159
1160     Float_t xstart = etp.GetX();
1161     if(xstart>3.) {
1162       d0z0[0]=-999.;
1163       d0z0[1]=-999.;
1164       //printf("This method can be used only for propagation inside the beam pipe \n");
1165       return kFALSE;
1166     }
1167
1168
1169     AliAODVertex *vtx =(AliAODVertex*)(evt->GetPrimaryVertex());
1170     Double_t fBzkG = evt->GetMagneticField(); // z componenent of field in kG
1171     ok = etp.PropagateToDCA(vtx,fBzkG,kVeryBig,d0z0,covd0z0);
1172     //ok = copy.PropagateToDCA(vtx,fBzkG,kVeryBig,d0z0,covd0z0);
1173   }
1174   if(!ok){
1175     d0z0[0]=-999.;
1176     d0z0[1]=-999.;
1177   }
1178   return ok;
1179 }
1180
1181
1182 Bool_t AlidNdPtAnalysisPbPbAOD::IsMCTrackAccepted(AliAODMCParticle *part)
1183 {
1184   if(!part) return kFALSE;
1185
1186   Double_t charge = part->Charge()/3.;
1187   if (TMath::Abs(charge) < 0.001) return kFALSE;
1188
1189   return kTRUE;
1190 }
1191
1192 const char * AlidNdPtAnalysisPbPbAOD::GetParticleName(Int_t pdg)
1193 {
1194   TParticlePDG * p1 = TDatabasePDG::Instance()->GetParticle(pdg);
1195   if(p1) return p1->GetName();
1196   return Form("%d", pdg);
1197 }
1198
1199 AliGenHijingEventHeader* AlidNdPtAnalysisPbPbAOD::GetHijingEventHeader(AliAODMCHeader *header)
1200 {
1201   //
1202   // inspired by PWGJE/AliPWG4HighPtSpectra.cxx
1203   //
1204
1205   if(!header) return 0x0;
1206   AliGenHijingEventHeader* hijingGenHeader = NULL;
1207
1208   TList* headerList = header->GetCocktailHeaders();
1209
1210   for(Int_t i = 0; i < headerList->GetEntries(); i++)
1211   {
1212     hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(headerList->At(i));
1213     if(hijingGenHeader) break;
1214   }
1215
1216   if(!hijingGenHeader) return 0x0;
1217
1218   return hijingGenHeader;
1219 }
1220
1221 AliGenPythiaEventHeader* AlidNdPtAnalysisPbPbAOD::GetPythiaEventHeader(AliAODMCHeader *header)
1222 {
1223   //
1224   // inspired by PWGJE/AliPWG4HighPtSpectra.cxx
1225   //
1226
1227   if(!header) return 0x0;
1228   AliGenPythiaEventHeader* PythiaGenHeader = NULL;
1229
1230   TList* headerList = header->GetCocktailHeaders();
1231
1232   for(Int_t i = 0; i < headerList->GetEntries(); i++)
1233   {
1234     PythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
1235     if(PythiaGenHeader) break;
1236   }
1237
1238   if(!PythiaGenHeader) return 0x0;
1239
1240   return PythiaGenHeader;
1241 }
1242
1243 //________________________________________________________________________
1244 Bool_t AlidNdPtAnalysisPbPbAOD::IsHijingParticle(const AliAODMCParticle *part, AliGenHijingEventHeader* hijingGenHeader){
1245
1246   // Check whether a particle is from Hijing or some injected
1247   // returns kFALSE if particle is injected
1248
1249   if(part->Label() > (hijingGenHeader->NProduced()-1)) return kFALSE;
1250   return kTRUE;
1251 }
1252
1253 //________________________________________________________________________
1254 Bool_t AlidNdPtAnalysisPbPbAOD::IsPythiaParticle(const AliAODMCParticle *part, AliGenPythiaEventHeader* pythiaGenHeader){
1255
1256   // Check whether a particle is from Pythia or some injected
1257
1258   if(part->Label() > (pythiaGenHeader->NProduced()-1)) return kFALSE;
1259   return kTRUE;
1260 }
1261
1262 Double_t* AlidNdPtAnalysisPbPbAOD::GetArrayClone(Int_t n, Double_t* source)
1263 {
1264   if (!source || n==0) return 0;
1265   Double_t* dest = new Double_t[n];
1266   for (Int_t i=0; i<n ; i++) { dest[i] = source[i]; }
1267   return dest;
1268 }
1269
1270 void AlidNdPtAnalysisPbPbAOD::Terminate(Option_t *)
1271 {
1272
1273 }
1274
1275