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