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 //------------------------------------------------------------------------------
  20
  21
  22 #include "AlidNdPtAnalysisPbPbAOD.h"
  23
  24
  25 using namespace std;
  26
  27 ClassImp(AlidNdPtAnalysisPbPbAOD)
  28
  29 // dummy constructor
  30 AlidNdPtAnalysisPbPbAOD::AlidNdPtAnalysisPbPbAOD() : AliAnalysisTaskSE(),
  31 fOutputList(0),
  32 // Histograms
  33 hPt(0),
  34 hMCPt(0),
  35 hnZvPtEtaCent(0),
  36 hnMCRecPrimZvPtEtaCent(0),
  37 hnMCGenZvPtEtaCent(0),
  38 hnMCRecSecZvPtEtaCent(0),
  39 hEventStatistics(0),
  40 hEventStatisticsCentrality(0),
  41 hMCEventStatisticsCentrality(0),
  42 hAllEventStatisticsCentrality(0),
  43 hEventStatisticsCentralityTrigger(0),
  44 hnZvMultCent(0),
  45 hTriggerStatistics(0),
  46 hMCTrackPdgCode(0),
  47 hMCTrackStatusCode(0),
  48 hCharge(0),
  49 hMCCharge(0),
  50 hMCPdgPt(0),
  51 hMCHijingPrim(0),
  52 hDCAPtAll(0),
  53 hDCAPtAccepted(0),
  54 hMCDCAPtSecondary(0),
  55 hMCDCAPtPrimary(0),
  56 hnCrossedRowsClustersChiPtEtaPhiAll(0),
  57 hnCrossedRowsClustersChiPtEtaPhiAcc(0),
  58 //global
  59 bIsMonteCarlo(0),
  60 // event cut variables
  61 dCutMaxZVertex(10.),
  62 // track kinematic cut variables
  63 dCutPtMin(0.15),
  64 dCutPtMax(1000.),
  65 dCutEtaMin(-0.8),
  66 dCutEtaMax(0.8),
  67 // track quality cut variables
  68 bCutRequireTPCRefit(kTRUE),
  69 dCutMinNumberOfCrossedRows(120.),
  70 dCutMinRatioCrossedRowsOverFindableClustersTPC(0.8),
  71 dCutMaxChi2PerClusterTPC(4.),
  72 dCutMaxFractionSharedTPCClusters(0.4),
  73 dCutMaxDCAToVertexZ(3.0),
  74 dCutMaxDCAToVertexXY(3.0),
  75 bCutRequireITSRefit(kTRUE),
  76 dCutMaxChi2PerClusterITS(36.),
  77 dCutDCAToVertex2D(kFALSE),
  78 dCutRequireSigmaToVertex(kFALSE),
  79 dCutMaxDCAToVertexXYPtDepPar0(0.0182),
  80 dCutMaxDCAToVertexXYPtDepPar1(0.0350),
  81 dCutMaxDCAToVertexXYPtDepPar2(1.01),
  82 bCutAcceptKinkDaughters(kFALSE),
  83 dCutMaxChi2TPCConstrainedGlobal(36.),
  84 // binning for THnSparse
  85 fMultNbins(0),
  86 fPtNbins(0),
  87 fPtCorrNbins(0),
  88 fEtaNbins(0),
  89 fZvNbins(0),
  90 fCentralityNbins(0),
  91 fBinsMult(0),
  92 fBinsPt(0),
  93 fBinsPtCorr(0),
  94 fBinsEta(0),
  95 fBinsZv(0),
  96 fBinsCentrality(0)
  97 {
  98
  99   fMultNbins = 0;
 100   fPtNbins = 0;
 101   fPtCorrNbins = 0;
 102   fEtaNbins = 0;
 103   fZvNbins = 0;
 104   fCentralityNbins = 0;
 105   fBinsMult = 0;
 106   fBinsPt = 0;
 107   fBinsPtCorr = 0;
 108   fBinsEta = 0;
 109   fBinsEta = 0;
 110   fBinsZv = 0;
 111   fBinsCentrality = 0;
 112
 113 }
 114
 115 AlidNdPtAnalysisPbPbAOD::AlidNdPtAnalysisPbPbAOD(const char *name) : AliAnalysisTaskSE(name),
 116 fOutputList(0),
 117 // Histograms
 118 hPt(0),
 119 hMCPt(0),
 120 hnZvPtEtaCent(0),
 121 hnMCRecPrimZvPtEtaCent(0),
 122 hnMCGenZvPtEtaCent(0),
 123 hnMCRecSecZvPtEtaCent(0),
 124 hEventStatistics(0),
 125 hEventStatisticsCentrality(0),
 126 hMCEventStatisticsCentrality(0),
 127 hAllEventStatisticsCentrality(0),
 128 hEventStatisticsCentralityTrigger(0),
 129 hnZvMultCent(0),
 130 hTriggerStatistics(0),
 131 hMCTrackPdgCode(0),
 132 hMCTrackStatusCode(0),
 133 hCharge(0),
 134 hMCCharge(0),
 135 hMCPdgPt(0),
 136 hMCHijingPrim(0),
 137 hDCAPtAll(0),
 138 hDCAPtAccepted(0),
 139 hMCDCAPtSecondary(0),
 140 hMCDCAPtPrimary(0),
 141 hnCrossedRowsClustersChiPtEtaPhiAll(0),
 142 hnCrossedRowsClustersChiPtEtaPhiAcc(0),
 143 //global
 144 bIsMonteCarlo(0),
 145 // event cut variables
 146 dCutMaxZVertex(10.),
 147 // track kinematic cut variables
 148 dCutPtMin(0.15),
 149 dCutPtMax(200.),
 150 dCutEtaMin(-0.8),
 151 dCutEtaMax(0.8),
 152 // track quality cut variables
 153 bCutRequireTPCRefit(kTRUE),
 154 dCutMinNumberOfCrossedRows(120.),
 155 dCutMinRatioCrossedRowsOverFindableClustersTPC(0.8),
 156 dCutMaxChi2PerClusterTPC(4.),
 157 dCutMaxFractionSharedTPCClusters(0.4),
 158 dCutMaxDCAToVertexZ(3.0),
 159 dCutMaxDCAToVertexXY(3.0),
 160 bCutRequireITSRefit(kTRUE),
 161 dCutMaxChi2PerClusterITS(36.),
 162 dCutDCAToVertex2D(kFALSE),
 163 dCutRequireSigmaToVertex(kFALSE),
 164 dCutMaxDCAToVertexXYPtDepPar0(0.0182),
 165 dCutMaxDCAToVertexXYPtDepPar1(0.0350),
 166 dCutMaxDCAToVertexXYPtDepPar2(1.01),
 167 bCutAcceptKinkDaughters(kFALSE),
 168 dCutMaxChi2TPCConstrainedGlobal(36.),
 169 // binning for THnSparse
 170 fMultNbins(0),
 171 fPtNbins(0),
 172 fPtCorrNbins(0),
 173 fEtaNbins(0),
 174 fZvNbins(0),
 175 fCentralityNbins(0),
 176 fBinsMult(0),
 177 fBinsPt(0),
 178 fBinsPtCorr(0),
 179 fBinsEta(0),
 180 fBinsZv(0),
 181 fBinsCentrality(0)
 182 {
 183   fMultNbins = 0;
 184   fPtNbins = 0;
 185   fPtCorrNbins = 0;
 186   fEtaNbins = 0;
 187   fZvNbins = 0;
 188   fCentralityNbins = 0;
 189   fBinsMult = 0;
 190   fBinsPt = 0;
 191   fBinsPtCorr = 0;
 192   fBinsEta = 0;
 193   fBinsEta = 0;
 194   fBinsZv = 0;
 195   fBinsCentrality = 0;
 196
 197   DefineOutput(1, TList::Class());
 198 }
 199
 200 // destructor
 201 AlidNdPtAnalysisPbPbAOD::~AlidNdPtAnalysisPbPbAOD()
 202 {
 203
 204   if(hnZvPtEtaCent) delete hnZvPtEtaCent; hnZvPtEtaCent = 0;
 205   if(hPt) delete hPt; hPt = 0;
 206   if(hnMCRecPrimZvPtEtaCent) delete hnMCRecPrimZvPtEtaCent; hnMCRecPrimZvPtEtaCent = 0;
 207   if(hnMCGenZvPtEtaCent) delete hnMCGenZvPtEtaCent; hnMCGenZvPtEtaCent = 0;
 208   if(hnMCRecSecZvPtEtaCent) delete hnMCRecSecZvPtEtaCent; hnMCRecSecZvPtEtaCent = 0;
 209   if(hMCPt) delete hMCPt; hMCPt = 0;
 210   if(hEventStatistics) delete hEventStatistics; hEventStatistics = 0;
 211   if(hEventStatisticsCentrality) delete hEventStatisticsCentrality; hEventStatisticsCentrality = 0;
 212   if(hMCEventStatisticsCentrality) delete hMCEventStatisticsCentrality; hMCEventStatisticsCentrality = 0;
 213   if(hAllEventStatisticsCentrality) delete hAllEventStatisticsCentrality; hAllEventStatisticsCentrality = 0;
 214   if(hEventStatisticsCentralityTrigger) delete hEventStatisticsCentralityTrigger; hEventStatisticsCentralityTrigger = 0;
 215   if(hnZvMultCent) delete hnZvMultCent; hnZvMultCent = 0;
 216   if(hTriggerStatistics) delete hTriggerStatistics; hTriggerStatistics = 0;
 217   if(hMCTrackPdgCode) delete hMCTrackPdgCode; hMCTrackPdgCode = 0;
 218   if(hMCTrackStatusCode) delete hMCTrackStatusCode; hMCTrackStatusCode = 0;
 219   if(hCharge) delete hCharge; hCharge = 0;
 220   if(hMCCharge) delete hMCCharge; hMCCharge = 0;
 221   if(hMCPdgPt) delete hMCPdgPt; hMCPdgPt = 0;
 222   if(hMCHijingPrim) delete hMCHijingPrim; hMCHijingPrim = 0;
 223   if(hDCAPtAll) delete hDCAPtAll; hDCAPtAll = 0;
 224   if(hDCAPtAccepted) delete hDCAPtAccepted; hDCAPtAccepted = 0;
 225   if(hMCDCAPtSecondary) delete hMCDCAPtSecondary; hMCDCAPtSecondary = 0;
 226   if(hMCDCAPtPrimary) delete hMCDCAPtPrimary; hMCDCAPtPrimary = 0;
 227   if(hMCDCAPtSecondary) delete hMCDCAPtSecondary; hMCDCAPtSecondary = 0;
 228   if(hMCDCAPtPrimary) delete hMCDCAPtPrimary; hMCDCAPtPrimary = 0;
 229   if(hnCrossedRowsClustersChiPtEtaPhiAll) delete hnCrossedRowsClustersChiPtEtaPhiAll; hnCrossedRowsClustersChiPtEtaPhiAll = 0;
 230   if(hnCrossedRowsClustersChiPtEtaPhiAcc) delete hnCrossedRowsClustersChiPtEtaPhiAcc; hnCrossedRowsClustersChiPtEtaPhiAcc = 0;
 231 }
 232
 233 void AlidNdPtAnalysisPbPbAOD::UserCreateOutputObjects()
 234 {
 235   // create all output histograms here
 236   OpenFile(1, "RECREATE");
 237
 238   fOutputList = new TList();
 239   fOutputList->SetOwner();
 240
 241   //define default binning
 242   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 };
 243   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};
 244   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};
 245   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};
 246   Double_t binsZvDefault[13] = {-30.,-25.,-20.,-15.,-10.,-5.,0.,5.,10.,15.,20.,25.,30.};
 247   Double_t binsCentralityDefault[12] = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.};
 248
 249   // if no binning is set, use the default
 250   if (!fBinsMult)       { SetBinsMult(48,binsMultDefault); }
 251   if (!fBinsPt)         { SetBinsPt(82,binsPtDefault); }
 252   if (!fBinsPtCorr)     { SetBinsPtCorr(37,binsPtCorrDefault); }
 253   if (!fBinsEta)        { SetBinsEta(31,binsEtaDefault); }
 254   if (!fBinsZv)         { SetBinsZv(13,binsZvDefault); }
 255   if (!fBinsCentrality) { SetBinsCentrality(12,binsCentralityDefault); }
 256
 257   Int_t binsZvPtEtaCent[4]={fZvNbins-1,fPtNbins-1,fEtaNbins-1,fCentralityNbins-1};
 258   Int_t binsZvMultCent[3]={fZvNbins-1,fMultNbins-1,fCentralityNbins-1};
 259
 260   // define Histograms
 261   hnZvPtEtaCent = new THnSparseF("hnZvPtEtaCent","Zv:Pt:Eta:Centrality",4,binsZvPtEtaCent);
 262   hnZvPtEtaCent->SetBinEdges(0,fBinsZv);
 263   hnZvPtEtaCent->SetBinEdges(1,fBinsPt);
 264   hnZvPtEtaCent->SetBinEdges(2,fBinsEta);
 265   hnZvPtEtaCent->SetBinEdges(3,fBinsCentrality);
 266   hnZvPtEtaCent->GetAxis(0)->SetTitle("Zv (cm)");
 267   hnZvPtEtaCent->GetAxis(1)->SetTitle("Pt (GeV/c)");
 268   hnZvPtEtaCent->GetAxis(2)->SetTitle("Eta");
 269   hnZvPtEtaCent->GetAxis(3)->SetTitle("Centrality");
 270   hnZvPtEtaCent->Sumw2();
 271
 272   hnMCRecPrimZvPtEtaCent = new THnSparseF("hnMCRecPrimZvPtEtaCent","mcZv:mcPt:mcEta:Centrality",4,binsZvPtEtaCent);
 273   hnMCRecPrimZvPtEtaCent->SetBinEdges(0,fBinsZv);
 274   hnMCRecPrimZvPtEtaCent->SetBinEdges(1,fBinsPt);
 275   hnMCRecPrimZvPtEtaCent->SetBinEdges(2,fBinsEta);
 276   hnMCRecPrimZvPtEtaCent->SetBinEdges(3,fBinsCentrality);
 277   hnMCRecPrimZvPtEtaCent->GetAxis(0)->SetTitle("MC Zv (cm)");
 278   hnMCRecPrimZvPtEtaCent->GetAxis(1)->SetTitle("MC Pt (GeV/c)");
 279   hnMCRecPrimZvPtEtaCent->GetAxis(2)->SetTitle("MC Eta");
 280   hnMCRecPrimZvPtEtaCent->GetAxis(3)->SetTitle("Centrality");
 281   hnMCRecPrimZvPtEtaCent->Sumw2();
 282
 283   hnMCGenZvPtEtaCent = new THnSparseF("hnMCGenZvPtEtaCent","mcZv:mcPt:mcEta:Centrality",4,binsZvPtEtaCent);
 284   hnMCGenZvPtEtaCent->SetBinEdges(0,fBinsZv);
 285   hnMCGenZvPtEtaCent->SetBinEdges(1,fBinsPt);
 286   hnMCGenZvPtEtaCent->SetBinEdges(2,fBinsEta);
 287   hnMCGenZvPtEtaCent->SetBinEdges(3,fBinsCentrality);
 288   hnMCGenZvPtEtaCent->GetAxis(0)->SetTitle("MC Zv (cm)");
 289   hnMCGenZvPtEtaCent->GetAxis(1)->SetTitle("MC Pt (GeV/c)");
 290   hnMCGenZvPtEtaCent->GetAxis(2)->SetTitle("MC Eta");
 291   hnMCGenZvPtEtaCent->GetAxis(3)->SetTitle("Centrality");
 292   hnMCGenZvPtEtaCent->Sumw2();
 293
 294   hnMCRecSecZvPtEtaCent = new THnSparseF("hnMCRecSecZvPtEtaCent","mcZv:mcPt:mcEta:Centrality",4,binsZvPtEtaCent);
 295   hnMCRecSecZvPtEtaCent->SetBinEdges(0,fBinsZv);
 296   hnMCRecSecZvPtEtaCent->SetBinEdges(1,fBinsPt);
 297   hnMCRecSecZvPtEtaCent->SetBinEdges(2,fBinsEta);
 298   hnMCRecSecZvPtEtaCent->SetBinEdges(3,fBinsCentrality);
 299   hnMCRecSecZvPtEtaCent->GetAxis(0)->SetTitle("MC Sec Zv (cm)");
 300   hnMCRecSecZvPtEtaCent->GetAxis(1)->SetTitle("MC Sec Pt (GeV/c)");
 301   hnMCRecSecZvPtEtaCent->GetAxis(2)->SetTitle("MC Sec Eta");
 302   hnMCRecSecZvPtEtaCent->GetAxis(3)->SetTitle("Centrality");
 303   hnMCRecSecZvPtEtaCent->Sumw2();
 304
 305   hPt = new TH1F("hPt","hPt",2000,0,200);
 306   hPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 307   hPt->GetYaxis()->SetTitle("dN/dp_{T}");
 308   hPt->Sumw2();
 309
 310   hMCPt = new TH1F("hMCPt","hMCPt",2000,0,200);
 311   hMCPt->GetXaxis()->SetTitle("MC p_{T} (GeV/c)");
 312   hMCPt->GetYaxis()->SetTitle("dN/dp_{T}");
 313   hMCPt->Sumw2();
 314
 315   hEventStatistics = new TH1F("hEventStatistics","hEventStatistics",10,0,10);
 316   hEventStatistics->GetYaxis()->SetTitle("number of events");
 317   hEventStatistics->SetBit(TH1::kCanRebin);
 318
 319   hEventStatisticsCentrality = new TH1F("hEventStatisticsCentrality","hEventStatisticsCentrality",fCentralityNbins-1, fBinsCentrality);
 320   hEventStatisticsCentrality->GetYaxis()->SetTitle("number of events");
 321
 322   hMCEventStatisticsCentrality = new TH1F("hMCEventStatisticsCentrality","hMCEventStatisticsCentrality",fCentralityNbins-1, fBinsCentrality);
 323   hMCEventStatisticsCentrality->GetYaxis()->SetTitle("number of MC events");
 324
 325   hAllEventStatisticsCentrality = new TH1F("hAllEventStatisticsCentrality","hAllEventStatisticsCentrality",fCentralityNbins-1, fBinsCentrality);
 326   hAllEventStatisticsCentrality->GetYaxis()->SetTitle("number of events");
 327
 328   hEventStatisticsCentralityTrigger = new TH2F("hEventStatisticsCentralityTrigger","hEventStatisticsCentralityTrigger;centrality;trigger",100,0,100,3,0,3);
 329   hEventStatisticsCentralityTrigger->Sumw2();
 330
 331   hnZvMultCent = new THnSparseF("hnZvMultCent","Zv:mult:Centrality",3,binsZvMultCent);
 332   hnZvMultCent->SetBinEdges(0,fBinsZv);
 333   hnZvMultCent->SetBinEdges(1,fBinsMult);
 334   hnZvMultCent->SetBinEdges(2,fBinsCentrality);
 335   hnZvMultCent->GetAxis(0)->SetTitle("Zv (cm)");
 336   hnZvMultCent->GetAxis(1)->SetTitle("N_{acc}");
 337   hnZvMultCent->GetAxis(2)->SetTitle("Centrality");
 338   hnZvMultCent->Sumw2();
 339
 340   hTriggerStatistics = new TH1F("hTriggerStatistics","hTriggerStatistics",10,0,10);
 341   hTriggerStatistics->GetYaxis()->SetTitle("number of events");
 342
 343   hMCTrackPdgCode = new TH1F("hMCTrackPdgCode","hMCTrackPdgCode",100,0,10);
 344   hMCTrackPdgCode->GetYaxis()->SetTitle("number of tracks");
 345   hMCTrackPdgCode->SetBit(TH1::kCanRebin);
 346
 347   hMCTrackStatusCode = new TH1F("hMCTrackStatusCode","hMCTrackStatusCode",100,0,10);
 348   hMCTrackStatusCode->GetYaxis()->SetTitle("number of tracks");
 349   hMCTrackStatusCode->SetBit(TH1::kCanRebin);
 350
 351   hCharge = new TH1F("hCharge","hCharge",30, -5, 5);
 352   hCharge->GetXaxis()->SetTitle("Charge");
 353   hCharge->GetYaxis()->SetTitle("number of tracks");
 354
 355   hMCCharge = new TH1F("hMCCharge","hMCCharge",30, -5, 5);
 356   hMCCharge->GetXaxis()->SetTitle("MC Charge");
 357   hMCCharge->GetYaxis()->SetTitle("number of tracks");
 358
 359   hMCPdgPt = new TH2F("hMCPdgPt","hMCPdgPt",fPtNbins-1, fBinsPt, 100,0,100);
 360   hMCPdgPt->GetYaxis()->SetTitle("particle");
 361   hMCPdgPt->GetXaxis()->SetTitle("Pt (GeV/c)");
 362
 363   hMCHijingPrim = new TH1F("hMCHijingPrim","hMCHijingPrim",2,0,2);
 364   hMCPdgPt->GetYaxis()->SetTitle("number of particles");
 365
 366
 367
 368   Int_t binsDCAxyDCAzPtEtaPhi[6] = { 200,200, fPtNbins-1, fEtaNbins-1, 36, fCentralityNbins-1};
 369   Double_t minDCAxyDCAzPtEtaPhi[6] = { -5, -5, 0, -1.5, 0., 0, };
 370   Double_t maxDCAxyDCAzPtEtaPhi[6] = { 5., 5., 100, 1.5, 2.*TMath::Pi(), 100};
 371
 372   hDCAPtAll = new THnSparseF("hDCAPtAll","hDCAPtAll",6, binsDCAxyDCAzPtEtaPhi, minDCAxyDCAzPtEtaPhi, maxDCAxyDCAzPtEtaPhi);
 373   hDCAPtAccepted = new THnSparseF("hDCAPtAccepted","hDCAPtAccepted",6, binsDCAxyDCAzPtEtaPhi, minDCAxyDCAzPtEtaPhi, maxDCAxyDCAzPtEtaPhi);
 374   hMCDCAPtSecondary = new THnSparseF("hMCDCAPtSecondary","hMCDCAPtSecondary",6, binsDCAxyDCAzPtEtaPhi, minDCAxyDCAzPtEtaPhi, maxDCAxyDCAzPtEtaPhi);
 375   hMCDCAPtPrimary = new THnSparseF("hMCDCAPtPrimary","hMCDCAPtPrimary",6, binsDCAxyDCAzPtEtaPhi, minDCAxyDCAzPtEtaPhi, maxDCAxyDCAzPtEtaPhi);
 376
 377   hDCAPtAll->SetBinEdges(2, fBinsPt);
 378   hDCAPtAccepted->SetBinEdges(2, fBinsPt);
 379   hMCDCAPtSecondary->SetBinEdges(2, fBinsPt);
 380   hMCDCAPtPrimary->SetBinEdges(2, fBinsPt);
 381
 382   hDCAPtAll->SetBinEdges(3, fBinsEta);
 383   hDCAPtAccepted->SetBinEdges(3, fBinsEta);
 384   hMCDCAPtSecondary->SetBinEdges(3, fBinsEta);
 385   hMCDCAPtPrimary->SetBinEdges(3, fBinsEta);
 386
 387   hDCAPtAll->SetBinEdges(5, fBinsCentrality);
 388   hDCAPtAccepted->SetBinEdges(5, fBinsCentrality);
 389   hMCDCAPtSecondary->SetBinEdges(5, fBinsCentrality);
 390   hMCDCAPtPrimary->SetBinEdges(5, fBinsCentrality);
 391
 392   hDCAPtAll->Sumw2();
 393   hDCAPtAccepted->Sumw2();
 394   hMCDCAPtSecondary->Sumw2();
 395   hMCDCAPtPrimary->Sumw2();
 396
 397   hDCAPtAll->GetAxis(0)->SetTitle("DCA_{xy} (cm)");
 398   hDCAPtAll->GetAxis(1)->SetTitle("DCA_{z} (cm)");
 399   hDCAPtAll->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
 400   hDCAPtAll->GetAxis(3)->SetTitle("#eta");
 401   hDCAPtAll->GetAxis(4)->SetTitle("#phi");
 402   hDCAPtAll->GetAxis(5)->SetTitle("Centrality");
 403
 404   hDCAPtAccepted->GetAxis(0)->SetTitle("DCA_{xy} (cm)");
 405   hDCAPtAccepted->GetAxis(1)->SetTitle("DCA_{z} (cm)");
 406   hDCAPtAccepted->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
 407   hDCAPtAccepted->GetAxis(3)->SetTitle("#eta");
 408   hDCAPtAccepted->GetAxis(4)->SetTitle("#phi");
 409   hDCAPtAccepted->GetAxis(5)->SetTitle("Centrality");
 410
 411   hMCDCAPtSecondary->GetAxis(0)->SetTitle("DCA_{xy} (cm)");
 412   hMCDCAPtSecondary->GetAxis(1)->SetTitle("DCA_{z} (cm)");
 413   hMCDCAPtSecondary->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
 414   hMCDCAPtSecondary->GetAxis(3)->SetTitle("#eta");
 415   hMCDCAPtSecondary->GetAxis(4)->SetTitle("#phi");
 416   hMCDCAPtSecondary->GetAxis(5)->SetTitle("Centrality");
 417
 418   hMCDCAPtPrimary->GetAxis(0)->SetTitle("DCA_{xy} (cm)");
 419   hMCDCAPtPrimary->GetAxis(1)->SetTitle("DCA_{z} (cm)");
 420   hMCDCAPtPrimary->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
 421   hMCDCAPtPrimary->GetAxis(3)->SetTitle("#eta");
 422   hMCDCAPtPrimary->GetAxis(4)->SetTitle("#phi");
 423   hMCDCAPtPrimary->GetAxis(5)->SetTitle("Centrality");
 424
 425   Int_t binsRowsClusterChiPtEtaPhi[7] = { 32,32, 100, fPtNbins-1, fEtaNbins-1, 36, fCentralityNbins-1};
 426   Double_t minRowsClusterChiPtEtaPhi[7] = { 0, 0, 0, 0, -1.5, 0., 0.};
 427   Double_t maxRowsClusterChiPtEtaPhi[7] = { 159, 159, 10, 100, 1.5, 2.*TMath::Pi(), 100.};
 428
 429   hnCrossedRowsClustersChiPtEtaPhiAll = new THnSparseF("hnCrossedRowsClustersChiPtEtaPhiAll","hnCrossedRowsClustersChiPtEtaPhiAll",7,binsRowsClusterChiPtEtaPhi, minRowsClusterChiPtEtaPhi, maxRowsClusterChiPtEtaPhi);
 430
 431   hnCrossedRowsClustersChiPtEtaPhiAcc = new THnSparseF("hnCrossedRowsClustersChiPtEtaPhiAcc","hnCrossedRowsClustersChiPtEtaPhiAcc",7,binsRowsClusterChiPtEtaPhi, minRowsClusterChiPtEtaPhi, maxRowsClusterChiPtEtaPhi);
 432
 433   hnCrossedRowsClustersChiPtEtaPhiAll->Sumw2();
 434   hnCrossedRowsClustersChiPtEtaPhiAll->SetBinEdges(3, fBinsPt);
 435   hnCrossedRowsClustersChiPtEtaPhiAll->SetBinEdges(4, fBinsEta);
 436   hnCrossedRowsClustersChiPtEtaPhiAll->SetBinEdges(6, fBinsCentrality);
 437   hnCrossedRowsClustersChiPtEtaPhiAll->GetAxis(0)->SetTitle("NcrossedRows before Cut");
 438   hnCrossedRowsClustersChiPtEtaPhiAll->GetAxis(1)->SetTitle("Nclusters before Cut");
 439   hnCrossedRowsClustersChiPtEtaPhiAll->GetAxis(2)->SetTitle("#Chi^{2}/cluster before Cut");
 440   hnCrossedRowsClustersChiPtEtaPhiAll->GetAxis(3)->SetTitle("p_{T} (GeV/c)");
 441   hnCrossedRowsClustersChiPtEtaPhiAll->GetAxis(4)->SetTitle("#eta");
 442   hnCrossedRowsClustersChiPtEtaPhiAll->GetAxis(5)->SetTitle("#phi");
 443   hnCrossedRowsClustersChiPtEtaPhiAll->GetAxis(6)->SetTitle("Centrality");
 444
 445   hnCrossedRowsClustersChiPtEtaPhiAcc->Sumw2();
 446   hnCrossedRowsClustersChiPtEtaPhiAcc->SetBinEdges(3, fBinsPt);
 447   hnCrossedRowsClustersChiPtEtaPhiAcc->SetBinEdges(4, fBinsEta);
 448   hnCrossedRowsClustersChiPtEtaPhiAcc->SetBinEdges(6, fBinsCentrality);
 449   hnCrossedRowsClustersChiPtEtaPhiAcc->GetAxis(0)->SetTitle("NcrossedRows after Cut");
 450   hnCrossedRowsClustersChiPtEtaPhiAcc->GetAxis(1)->SetTitle("Nclusters after Cut");
 451   hnCrossedRowsClustersChiPtEtaPhiAcc->GetAxis(2)->SetTitle("#Chi^{2}/cluster after Cut");
 452   hnCrossedRowsClustersChiPtEtaPhiAcc->GetAxis(3)->SetTitle("p_{T} (GeV/c)");
 453   hnCrossedRowsClustersChiPtEtaPhiAcc->GetAxis(4)->SetTitle("#eta");
 454   hnCrossedRowsClustersChiPtEtaPhiAcc->GetAxis(5)->SetTitle("#phi");
 455   hnCrossedRowsClustersChiPtEtaPhiAcc->GetAxis(6)->SetTitle("Centrality");
 456
 457   // Add Histos, Profiles etc to List
 458   fOutputList->Add(hnZvPtEtaCent);
 459   fOutputList->Add(hPt);
 460   fOutputList->Add(hnMCRecPrimZvPtEtaCent);
 461   fOutputList->Add(hnMCGenZvPtEtaCent);
 462   fOutputList->Add(hnMCRecSecZvPtEtaCent);
 463   fOutputList->Add(hMCPt);
 464   fOutputList->Add(hEventStatistics);
 465   fOutputList->Add(hEventStatisticsCentrality);
 466   fOutputList->Add(hMCEventStatisticsCentrality);
 467   fOutputList->Add(hAllEventStatisticsCentrality);
 468   fOutputList->Add(hEventStatisticsCentralityTrigger);
 469   fOutputList->Add(hnZvMultCent);
 470   fOutputList->Add(hTriggerStatistics);
 471   fOutputList->Add(hMCTrackPdgCode);
 472   fOutputList->Add(hMCTrackStatusCode);
 473   fOutputList->Add(hCharge);
 474   fOutputList->Add(hMCCharge);
 475   fOutputList->Add(hMCPdgPt);
 476   fOutputList->Add(hMCHijingPrim);
 477   fOutputList->Add(hDCAPtAll);
 478   fOutputList->Add(hDCAPtAccepted);
 479   fOutputList->Add(hMCDCAPtSecondary);
 480   fOutputList->Add(hMCDCAPtPrimary);
 481   fOutputList->Add(hnCrossedRowsClustersChiPtEtaPhiAll);
 482   fOutputList->Add(hnCrossedRowsClustersChiPtEtaPhiAcc);
 483
 484   PostData(1, fOutputList);
 485 }
 486
 487 void AlidNdPtAnalysisPbPbAOD::UserExec(Option_t *option)
 488 {
 489
 490   // Main Loop
 491   // called for each event
 492   hEventStatistics->Fill("all events",1);
 493
 494   // set ZERO pointers:
 495   AliInputEventHandler *inputHandler = NULL;
 496   AliAODTrack *track = NULL;
 497   AliAODMCParticle *mcPart = NULL;
 498   AliAODMCHeader *mcHdr = NULL;
 499   AliGenHijingEventHeader *genHijingHeader = NULL;
 500   //AliGenPythiaEventHeader *genPythiaHeader = NULL;
 501
 502   Bool_t bIsEventSelectedMB = kFALSE;
 503   Bool_t bIsEventSelectedSemi = kFALSE;
 504   Bool_t bIsEventSelectedCentral = kFALSE;
 505   Bool_t bIsEventSelected = kFALSE;
 506   Bool_t bIsPrimary = kFALSE;
 507   Bool_t bIsHijingParticle = kFALSE;
 508   Bool_t bMotherIsHijingParticle = kFALSE;
 509   //Bool_t bIsPythiaParticle = kFALSE;
 510   Bool_t bEventHasATrack = kFALSE;
 511   Bool_t bEventHasATrackInRange = kFALSE;
 512   Int_t nTriggerFired = 0;
 513
 514
 515   Double_t dMCTrackZvPtEtaCent[4] = {0};
 516   Double_t dTrackZvPtEtaCent[4] = {0};
 517
 518   Double_t dDCA[2] = {0};
 519
 520   Double_t dMCEventZv = -100;
 521   Double_t dEventZv = -100;
 522   Int_t iAcceptedMultiplicity = 0;
 523
 524   bIsMonteCarlo = kFALSE;
 525
 526   AliAODEvent *eventAOD = 0x0;
 527   eventAOD = dynamic_cast<AliAODEvent*>( InputEvent() );
 528   if (!eventAOD) {
 529     AliWarning("ERROR: eventAOD not available \n");
 530     return;
 531   }
 532
 533   // check, which trigger has been fired
 534   inputHandler = (AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 535   bIsEventSelectedMB = ( inputHandler->IsEventSelected() & AliVEvent::kMB);
 536   bIsEventSelectedSemi = ( inputHandler->IsEventSelected() & AliVEvent::kSemiCentral);
 537   bIsEventSelectedCentral = ( inputHandler->IsEventSelected() & AliVEvent::kCentral);
 538
 539   if(bIsEventSelectedMB || bIsEventSelectedSemi || bIsEventSelectedCentral) hTriggerStatistics->Fill("all triggered events",1);
 540   if(bIsEventSelectedMB) { hTriggerStatistics->Fill("MB trigger",1); nTriggerFired++; }
 541   if(bIsEventSelectedSemi) { hTriggerStatistics->Fill("SemiCentral trigger",1); nTriggerFired++; }
 542   if(bIsEventSelectedCentral) { hTriggerStatistics->Fill("Central trigger",1); nTriggerFired++; }
 543   if(nTriggerFired == 0) { hTriggerStatistics->Fill("No trigger",1); }
 544
 545   bIsEventSelected = ( inputHandler->IsEventSelected() & GetCollisionCandidates() );
 546
 547   // only take tracks of events, which are triggered
 548   if(nTriggerFired == 0) { return; }
 549
 550   //   if( !bIsEventSelected || nTriggerFired>1 ) return;
 551
 552   //   hEventStatistics->Fill("events with only coll. cand.", 1);
 553
 554
 555
 556   // check if there is a stack, if yes, then do MC loop
 557   TList *list = eventAOD->GetList();
 558   TClonesArray *stack = 0x0;
 559   stack = (TClonesArray*)list->FindObject(AliAODMCParticle::StdBranchName());
 560
 561   if( stack )
 562   {
 563     bIsMonteCarlo = kTRUE;
 564
 565     mcHdr = (AliAODMCHeader*)list->FindObject(AliAODMCHeader::StdBranchName());
 566
 567     genHijingHeader = GetHijingEventHeader(mcHdr);
 568     //     genPythiaHeader = GetPythiaEventHeader(mcHdr);
 569
 570     if(!genHijingHeader) { return; }
 571
 572     //     if(!genPythiaHeader)  { return; }
 573
 574     dMCEventZv = mcHdr->GetVtxZ();
 575     dMCTrackZvPtEtaCent[0] = dMCEventZv;
 576     hEventStatistics->Fill("MC all events",1);
 577   }
 578
 579   AliCentrality* aCentrality = eventAOD->GetCentrality();
 580   Double_t dCentrality = aCentrality->GetCentralityPercentile("V0M");
 581
 582   if( dCentrality < 0 ) return;
 583   hEventStatistics->Fill("after centrality selection",1);
 584
 585
 586
 587   // start with MC truth analysis
 588   if(bIsMonteCarlo)
 589   {
 590
 591     if( dMCEventZv > dCutMaxZVertex )  { return; }
 592
 593     dMCTrackZvPtEtaCent[0] = dMCEventZv;
 594
 595     hEventStatistics->Fill("MC afterZv cut",1);
 596
 597     for(Int_t iMCtrack = 0; iMCtrack < stack->GetEntriesFast(); iMCtrack++)
 598     {
 599       mcPart =(AliAODMCParticle*)stack->At(iMCtrack);
 600
 601       // check for charge
 602       if( !(IsMCTrackAccepted(mcPart)) ) continue;
 603
 604       if(!IsHijingParticle(mcPart, genHijingHeader)) { continue; }
 605
 606       if(mcPart->IsPhysicalPrimary() )
 607       {
 608         hMCHijingPrim->Fill("IsPhysicalPrimary",1);
 609       }
 610       else
 611       {
 612         hMCHijingPrim->Fill("NOT a primary",1);
 613         continue;
 614       }
 615
 616
 617       //
 618       // ======================== fill histograms ========================
 619       dMCTrackZvPtEtaCent[1] = mcPart->Pt();
 620       dMCTrackZvPtEtaCent[2] = mcPart->Eta();
 621       dMCTrackZvPtEtaCent[3] = dCentrality;
 622
 623       bEventHasATrack = kTRUE;
 624
 625       hnMCGenZvPtEtaCent->Fill(dMCTrackZvPtEtaCent);
 626
 627       if( (dMCTrackZvPtEtaCent[1] > dCutPtMin) &&
 628         (dMCTrackZvPtEtaCent[1] < dCutPtMax) &&
 629         (dMCTrackZvPtEtaCent[2] > dCutEtaMin) &&
 630         (dMCTrackZvPtEtaCent[2] < dCutEtaMax) )
 631       {
 632         hMCPt->Fill(mcPart->Pt());
 633         hMCCharge->Fill(mcPart->Charge()/3.);
 634         bEventHasATrackInRange = kTRUE;
 635       }
 636
 637     }
 638   } // isMonteCarlo
 639   if(bEventHasATrack) { hEventStatistics->Fill("MC events with tracks",1); }
 640   if(bEventHasATrackInRange)
 641   {
 642     hEventStatistics->Fill("MC events with tracks in range",1);
 643     hMCEventStatisticsCentrality->Fill(dCentrality);
 644   }
 645   bEventHasATrack = kFALSE;
 646   bEventHasATrackInRange = kFALSE;
 647
 648
 649
 650   // Loop over recontructed tracks
 651
 652   dEventZv = eventAOD->GetPrimaryVertex()->GetZ();
 653   if( TMath::Abs(dEventZv) > dCutMaxZVertex ) return;
 654
 655   hAllEventStatisticsCentrality->Fill(dCentrality/*, nTriggerFired*/);
 656
 657   hEventStatistics->Fill("after Zv cut",1);
 658
 659   dTrackZvPtEtaCent[0] = dEventZv;
 660
 661   for(Int_t itrack = 0; itrack < eventAOD->GetNumberOfTracks(); itrack++)
 662   {
 663     track = eventAOD->GetTrack(itrack);
 664     if(!track) continue;
 665
 666     mcPart = NULL;
 667     dMCTrackZvPtEtaCent[1] = 0;
 668     dMCTrackZvPtEtaCent[2] = 0;
 669     dMCTrackZvPtEtaCent[3] = 0;
 670
 671     bIsPrimary = kFALSE;
 672
 673     GetDCA(track, eventAOD, dDCA);
 674
 675     Double_t dDCAxyDCAzPt[5] = { dDCA[0], dDCA[1], track->Pt(), track->Eta(), track->Phi() };
 676
 677     hDCAPtAll->Fill(dDCAxyDCAzPt);
 678
 679     if( !(IsTrackAccepted(track)) ) continue;
 680
 681     dTrackZvPtEtaCent[1] = track->Pt();
 682     dTrackZvPtEtaCent[2] = track->Eta();
 683     dTrackZvPtEtaCent[3] = dCentrality;
 684
 685     if( bIsMonteCarlo )
 686     {
 687       mcPart = (AliAODMCParticle*)stack->At(TMath::Abs(track->GetLabel()));
 688       if( !mcPart ) { continue; }
 689
 690       // check for charge
 691       if( !(IsMCTrackAccepted(mcPart)) ) {  continue; }
 692
 693       bIsHijingParticle = IsHijingParticle(mcPart, genHijingHeader);
 694       //       bIsPythiaParticle = IsPythiaParticle(mcPart, genPythiaHeader);
 695
 696       //       if(!bIsHijingParticle) continue; // only take real tracks, not injected ones
 697
 698       bIsPrimary = mcPart->IsPhysicalPrimary();
 699
 700       dMCTrackZvPtEtaCent[1] = mcPart->Pt();
 701       dMCTrackZvPtEtaCent[2] = mcPart->Eta();
 702       dMCTrackZvPtEtaCent[3] = dCentrality;
 703
 704       if(bIsPrimary && bIsHijingParticle)
 705       {
 706         hnMCRecPrimZvPtEtaCent->Fill(dMCTrackZvPtEtaCent);
 707         hMCDCAPtPrimary->Fill(dDCAxyDCAzPt);
 708       }
 709
 710       if(!bIsPrimary /*&& !bIsHijingParticle*/)
 711       {
 712         Int_t indexMoth = mcPart->GetMother();
 713         if(indexMoth >= 0)
 714         {
 715           AliAODMCParticle* moth = (AliAODMCParticle*)stack->At(indexMoth);
 716           bMotherIsHijingParticle = IsHijingParticle(moth, genHijingHeader);
 717
 718           if(bMotherIsHijingParticle) // only store secondaries, which come from a not embedded signal!
 719           {
 720             hMCTrackStatusCode->Fill(Form("%d",mcPart->GetStatus()), 1);
 721             if(TMath::Abs(mcPart->Eta()) < 0.8) { hMCPdgPt->Fill(mcPart->Pt(), Form("%s",GetParticleName(mcPart->GetPdgCode())), 1); }
 722
 723             hnMCRecSecZvPtEtaCent->Fill(dMCTrackZvPtEtaCent);
 724             hMCDCAPtSecondary->Fill(dDCAxyDCAzPt);
 725             hMCTrackPdgCode->Fill(Form("%s_H%i_H%i",GetParticleName(moth->GetPdgCode()),bMotherIsHijingParticle, bIsHijingParticle), 1);
 726             //    delete moth;
 727           }
 728         }
 729       }
 730     } // end isMonteCarlo
 731
 732     // ======================== fill histograms ========================
 733
 734     // only keep prim and sec from not embedded signal
 735     Bool_t bKeepMCTrack = kFALSE;
 736     if(bIsMonteCarlo)
 737     {
 738       if( (bIsHijingParticle && bIsPrimary) ^ (bMotherIsHijingParticle && !bIsPrimary) )
 739       {
 740         bKeepMCTrack = kTRUE;
 741       }
 742       else
 743       {
 744         continue;
 745       }
 746     }
 747
 748     bEventHasATrack = kTRUE;
 749
 750     hnZvPtEtaCent->Fill(dTrackZvPtEtaCent);
 751     hDCAPtAccepted->Fill(dDCAxyDCAzPt);
 752
 753     if( (dTrackZvPtEtaCent[1] > dCutPtMin) &&
 754       (dTrackZvPtEtaCent[1] < dCutPtMax) &&
 755       (dTrackZvPtEtaCent[2] > dCutEtaMin) &&
 756       (dTrackZvPtEtaCent[2] < dCutEtaMax) )
 757     {
 758       iAcceptedMultiplicity++;
 759       bEventHasATrackInRange = kTRUE;
 760       hPt->Fill(track->Pt());
 761       hCharge->Fill(track->Charge());
 762     }
 763   } // end track loop
 764
 765   if(bEventHasATrack) { hEventStatistics->Fill("events with tracks",1); bEventHasATrack = kFALSE; }
 766
 767   if(bEventHasATrackInRange)
 768   {
 769     hEventStatistics->Fill("events with tracks in range",1);
 770     hEventStatisticsCentrality->Fill(dCentrality);
 771     bEventHasATrackInRange = kFALSE;
 772
 773     if(bIsEventSelectedMB) hEventStatisticsCentralityTrigger->Fill(dCentrality, 0);
 774     if(bIsEventSelectedSemi) hEventStatisticsCentralityTrigger->Fill(dCentrality, 1);
 775     if(bIsEventSelectedCentral) hEventStatisticsCentralityTrigger->Fill(dCentrality, 2);
 776   }
 777
 778   Double_t dEventZvMultCent[3] = {dEventZv, iAcceptedMultiplicity, dCentrality};
 779   hnZvMultCent->Fill(dEventZvMultCent);
 780
 781   PostData(1, fOutputList);
 782
 783 }
 784
 785 Bool_t AlidNdPtAnalysisPbPbAOD::IsTrackAccepted(AliAODTrack *tr)
 786 {
 787   if(!tr) return kFALSE;
 788
 789   if(tr->Charge()==0) { return kFALSE; }
 790
 791   Double_t dNClustersTPC = tr->GetTPCNcls();
 792   Double_t dCrossedRowsTPC = tr->GetTPCClusterInfo(2,1);
 793   Double_t dChi2PerClusterTPC = (dNClustersTPC>0)?tr->Chi2perNDF()*(dNClustersTPC-5)/dNClustersTPC:-1.; // see AliDielectronVarManager.h
 794
 795   //   hAllCrossedRowsTPC->Fill(dCrossedRowsTPC);
 796
 797   Double_t dRowClusterChiPtEtaPhi[6] = {dCrossedRowsTPC, dNClustersTPC, dChi2PerClusterTPC, tr->Pt(), tr->Eta(), tr->Phi() };
 798   hnCrossedRowsClustersChiPtEtaPhiAll->Fill(dRowClusterChiPtEtaPhi);
 799
 800   // filter bit 5
 801   if(!(tr->TestFilterBit(AliAODTrack::kTrkGlobal)) ) { return kFALSE; }
 802
 803   // filter bit 4
 804   //   if(!(tr->TestFilterBit(AliAODTrack::kTrkGlobalNoDCA)) ) { return kFALSE; }
 805
 806   //   hFilterCrossedRowsTPC->Fill(dCrossedRowsTPC);
 807
 808   if(dCrossedRowsTPC < GetCutMinNCrossedRowsTPC()) { return kFALSE; }
 809
 810   hnCrossedRowsClustersChiPtEtaPhiAcc->Fill(dRowClusterChiPtEtaPhi);
 811
 812   //   hAccNclsTPC->Fill(dNClustersTPC);
 813   //   hAccCrossedRowsTPC->Fill(dCrossedRowsTPC);
 814   //   Double_t dFindableClustersTPC = tr->GetTPCNclsF();
 815   //   Double_t dChi2PerClusterTPC = (dNClustersTPC>0)?tr->Chi2perNDF()*(dNClustersTPC-5)/dNClustersTPC:-1.; // see AliDielectronVarManager.h
 816   //
 817   //   Bool_t bIsFromKink = kFALSE;
 818   //   if(tr->GetProdVertex()->GetType() == AliAODVertex::kKink) bIsFromKink = kTRUE;
 819   //
 820   //   // from AliAnalysisTaskPIDqa.cxx
 821   //   ULong_t uStatus = tr->GetStatus();
 822   //   Bool_t bHasRefitTPC = kFALSE;
 823   //   Bool_t bHasRefitITS = kFALSE;
 824   //
 825   //   if ((uStatus & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) bHasRefitTPC = kTRUE;
 826   //   if ((uStatus & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) bHasRefitITS = kTRUE;
 827   //
 828   //   // from AliDielectronVarManager.h
 829   //   Bool_t bHasHitInSPD = kFALSE;
 830   //   for (Int_t iC=0; iC<2; iC++)
 831   //   {
 832     //     if (((tr->GetITSClusterMap()) & (1<<(iC))) > 0) {  bHasHitInSPD = kTRUE;  }
 833     //   }
 834     //
 835     //   Double_t dNClustersITS = tr->GetITSNcls();
 836
 837     // cuts to be done:
 838     // TPC
 839     //   esdTrackCuts->SetMinNCrossedRowsTPC(70);
 840     //   esdTrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
 841     //
 842     //   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
 843     //   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 844     //   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 845     // ITS
 846     //   esdTrackCuts->SetRequireITSRefit(kTRUE);
 847     //   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,   AliESDtrackCuts::kAny);
 848     //
 849     //   esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0105+0.0350/pt^1.1");
 850     //   esdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
 851     //
 852     //   esdTrackCuts->SetMaxDCAToVertexZ(2);
 853     //   esdTrackCuts->SetDCAToVertex2D(kFALSE);
 854     //   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 855     //
 856     //   esdTrackCuts->SetMaxChi2PerClusterITS(36);
 857
 858
 859     return kTRUE;
 860 }
 861
 862 Bool_t AlidNdPtAnalysisPbPbAOD::GetDCA(const AliAODTrack *track, AliAODEvent *evt, Double_t d0z0[2])
 863 {
 864   // function adapted from AliDielectronVarManager.h
 865
 866   if(track->TestBit(AliAODTrack::kIsDCA)){
 867     d0z0[0]=track->DCA();
 868     d0z0[1]=track->ZAtDCA();
 869     return kTRUE;
 870   }
 871
 872   Bool_t ok=kFALSE;
 873   if(evt) {
 874     Double_t covd0z0[3];
 875     //AliAODTrack copy(*track);
 876     AliExternalTrackParam etp; etp.CopyFromVTrack(track);
 877
 878     Float_t xstart = etp.GetX();
 879     if(xstart>3.) {
 880       d0z0[0]=-999.;
 881       d0z0[1]=-999.;
 882       //printf("This method can be used only for propagation inside the beam pipe \n");
 883       return kFALSE;
 884     }
 885
 886
 887     AliAODVertex *vtx =(AliAODVertex*)(evt->GetPrimaryVertex());
 888     Double_t fBzkG = evt->GetMagneticField(); // z componenent of field in kG
 889     ok = etp.PropagateToDCA(vtx,fBzkG,kVeryBig,d0z0,covd0z0);
 890     //ok = copy.PropagateToDCA(vtx,fBzkG,kVeryBig,d0z0,covd0z0);
 891   }
 892   if(!ok){
 893     d0z0[0]=-999.;
 894     d0z0[1]=-999.;
 895   }
 896   return ok;
 897 }
 898
 899
 900 Bool_t AlidNdPtAnalysisPbPbAOD::IsMCTrackAccepted(AliAODMCParticle *part)
 901 {
 902   if(!part) return kFALSE;
 903
 904   Double_t charge = part->Charge()/3.;
 905   if (TMath::Abs(charge) < 0.001) return kFALSE;
 906
 907   return kTRUE;
 908 }
 909
 910 const char * AlidNdPtAnalysisPbPbAOD::GetParticleName(Int_t pdg)
 911 {
 912   TParticlePDG * p1 = TDatabasePDG::Instance()->GetParticle(pdg);
 913   if(p1) return p1->GetName();
 914   return Form("%d", pdg);
 915 }
 916
 917 AliGenHijingEventHeader* AlidNdPtAnalysisPbPbAOD::GetHijingEventHeader(AliAODMCHeader *header)
 918 {
 919   //
 920   // inspired by PWGJE/AliPWG4HighPtSpectra.cxx
 921   //
 922
 923   if(!header) return 0x0;
 924   AliGenHijingEventHeader* hijingGenHeader = NULL;
 925
 926   TList* headerList = header->GetCocktailHeaders();
 927
 928   for(Int_t i = 0; i < headerList->GetEntries(); i++)
 929   {
 930     hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(headerList->At(i));
 931     if(hijingGenHeader) break;
 932   }
 933
 934   if(!hijingGenHeader) return 0x0;
 935
 936   return hijingGenHeader;
 937 }
 938
 939 AliGenPythiaEventHeader* AlidNdPtAnalysisPbPbAOD::GetPythiaEventHeader(AliAODMCHeader *header)
 940 {
 941   //
 942   // inspired by PWGJE/AliPWG4HighPtSpectra.cxx
 943   //
 944
 945   if(!header) return 0x0;
 946   AliGenPythiaEventHeader* PythiaGenHeader = NULL;
 947
 948   TList* headerList = header->GetCocktailHeaders();
 949
 950   for(Int_t i = 0; i < headerList->GetEntries(); i++)
 951   {
 952     PythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
 953     if(PythiaGenHeader) break;
 954   }
 955
 956   if(!PythiaGenHeader) return 0x0;
 957
 958   return PythiaGenHeader;
 959 }
 960
 961 //________________________________________________________________________
 962 Bool_t AlidNdPtAnalysisPbPbAOD::IsHijingParticle(const AliAODMCParticle *part, AliGenHijingEventHeader* hijingGenHeader){
 963
 964   // Check whether a particle is from Hijing or some injected
 965
 966   if(part->Label() > (hijingGenHeader->NProduced()-1)) return kFALSE;
 967   return kTRUE;
 968 }
 969
 970 //________________________________________________________________________
 971 Bool_t AlidNdPtAnalysisPbPbAOD::IsPythiaParticle(const AliAODMCParticle *part, AliGenPythiaEventHeader* pythiaGenHeader){
 972
 973   // Check whether a particle is from Pythia or some injected
 974
 975   if(part->Label() > (pythiaGenHeader->NProduced()-1)) return kFALSE;
 976   return kTRUE;
 977 }
 978
 979 Double_t* AlidNdPtAnalysisPbPbAOD::GetArrayClone(Int_t n, Double_t* source)
 980 {
 981   if (!source || n==0) return 0;
 982   Double_t* dest = new Double_t[n];
 983   for (Int_t i=0; i<n ; i++) { dest[i] = source[i]; }
 984   return dest;
 985 }
 986
 987 void AlidNdPtAnalysisPbPbAOD::Terminate(Option_t *)
 988 {
 989
 990 }
 991
 992