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