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