PWGLF/SPECTRA/ChargedHadrons/dNdPt/AlidNdPtAnalysisPbPbAOD.cxx

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