PWGLF/SPECTRA/ChargedHadrons/dNdPt/AlidNdPtAnalysisPbPbAOD.cxx

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