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