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