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