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 //------------------------------------------------------------------------------
  20
  21
  22 #include "AlidNdPtAnalysisPbPbAOD.h"
  23
  24
  25 using namespace std;
  26
  27 ClassImp(AlidNdPtAnalysisPbPbAOD)
  28
  29 // dummy constructor
  30 AlidNdPtAnalysisPbPbAOD::AlidNdPtAnalysisPbPbAOD() : AliAnalysisTaskSE(),
  31 fOutputList(0),
  32 // Histograms
  33 hPt(0),
  34 hMCPt(0),
  35 hnZvPtEtaCent(0),
  36 hnMCRecPrimZvPtEtaCent(0),
  37 hnMCGenZvPtEtaCent(0),
  38 hnMCRecSecZvPtEtaCent(0),
  39 hEventStatistics(0),
  40 hEventStatisticsCentrality(0),
  41 hAllEventStatisticsCentrality(0),
  42 hnZvMultCent(0),
  43 hTriggerStatistics(0),
  44 hMCTrackPdgCode(0),
  45 hMCTrackStatusCode(0),
  46 hCharge(0),
  47 hMCCharge(0),
  48 hMCPdgPt(0),
  49 hMCHijingPrim(0),
  50 hAccNclsTPC(0),
  51 hAccCrossedRowsTPC(0),
  52 hDCAPtAll(0),
  53 hDCAPtAccepted(0),
  54 hMCDCAPtSecondary(0),
  55 hMCDCAPtPrimary(0),
  56 //global
  57 bIsMonteCarlo(0),
  58 // event cut variables
  59 dCutMaxZVertex(10.),
  60 // track kinematic cut variables
  61 dCutPtMin(0.15),
  62 dCutPtMax(1000.),
  63 dCutEtaMin(-0.8),
  64 dCutEtaMax(0.8),
  65 // track quality cut variables
  66 bCutRequireTPCRefit(kTRUE),
  67 dCutMinNumberOfCrossedRows(120.),
  68 dCutMinRatioCrossedRowsOverFindableClustersTPC(0.8),
  69 dCutMaxChi2PerClusterTPC(4.),
  70 dCutMaxFractionSharedTPCClusters(0.4),
  71 dCutMaxDCAToVertexZ(3.0),
  72 dCutMaxDCAToVertexXY(3.0),
  73 bCutRequireITSRefit(kTRUE),
  74 dCutMaxChi2PerClusterITS(36.),
  75 dCutDCAToVertex2D(kFALSE),
  76 dCutRequireSigmaToVertex(kFALSE),
  77 dCutMaxDCAToVertexXYPtDepPar0(0.0182),
  78 dCutMaxDCAToVertexXYPtDepPar1(0.0350),
  79 dCutMaxDCAToVertexXYPtDepPar2(1.01),
  80 bCutAcceptKinkDaughters(kFALSE),
  81 dCutMaxChi2TPCConstrainedGlobal(36.),
  82 // binning for THnSparse
  83 fMultNbins(0),
  84 fPtNbins(0),
  85 fPtCorrNbins(0),
  86 fEtaNbins(0),
  87 fZvNbins(0),
  88 fCentralityNbins(0),
  89 fBinsMult(0),
  90 fBinsPt(0),
  91 fBinsPtCorr(0),
  92 fBinsEta(0),
  93 fBinsZv(0),
  94 fBinsCentrality(0)
  95 {
  96
  97   fMultNbins = 0;
  98   fPtNbins = 0;
  99   fPtCorrNbins = 0;
 100   fEtaNbins = 0;
 101   fZvNbins = 0;
 102   fCentralityNbins = 0;
 103   fBinsMult = 0;
 104   fBinsPt = 0;
 105   fBinsPtCorr = 0;
 106   fBinsEta = 0;
 107   fBinsEta = 0;
 108   fBinsZv = 0;
 109   fBinsCentrality = 0;
 110
 111 }
 112
 113 AlidNdPtAnalysisPbPbAOD::AlidNdPtAnalysisPbPbAOD(const char *name) : AliAnalysisTaskSE(name),
 114 fOutputList(0),
 115 // Histograms
 116 hPt(0),
 117 hMCPt(0),
 118 hnZvPtEtaCent(0),
 119 hnMCRecPrimZvPtEtaCent(0),
 120 hnMCGenZvPtEtaCent(0),
 121 hnMCRecSecZvPtEtaCent(0),
 122 hEventStatistics(0),
 123 hEventStatisticsCentrality(0),
 124 hAllEventStatisticsCentrality(0),
 125 hnZvMultCent(0),
 126 hTriggerStatistics(0),
 127 hMCTrackPdgCode(0),
 128 hMCTrackStatusCode(0),
 129 hCharge(0),
 130 hMCCharge(0),
 131 hMCPdgPt(0),
 132 hMCHijingPrim(0),
 133 hAccNclsTPC(0),
 134 hAccCrossedRowsTPC(0),
 135 hDCAPtAll(0),
 136 hDCAPtAccepted(0),
 137 hMCDCAPtSecondary(0),
 138 hMCDCAPtPrimary(0),
 139 //global
 140 bIsMonteCarlo(0),
 141 // event cut variables
 142 dCutMaxZVertex(10.),
 143 // track kinematic cut variables
 144 dCutPtMin(0.15),
 145 dCutPtMax(200.),
 146 dCutEtaMin(-0.8),
 147 dCutEtaMax(0.8),
 148 // track quality cut variables
 149 bCutRequireTPCRefit(kTRUE),
 150 dCutMinNumberOfCrossedRows(120.),
 151 dCutMinRatioCrossedRowsOverFindableClustersTPC(0.8),
 152 dCutMaxChi2PerClusterTPC(4.),
 153 dCutMaxFractionSharedTPCClusters(0.4),
 154 dCutMaxDCAToVertexZ(3.0),
 155 dCutMaxDCAToVertexXY(3.0),
 156 bCutRequireITSRefit(kTRUE),
 157 dCutMaxChi2PerClusterITS(36.),
 158 dCutDCAToVertex2D(kFALSE),
 159 dCutRequireSigmaToVertex(kFALSE),
 160 dCutMaxDCAToVertexXYPtDepPar0(0.0182),
 161 dCutMaxDCAToVertexXYPtDepPar1(0.0350),
 162 dCutMaxDCAToVertexXYPtDepPar2(1.01),
 163 bCutAcceptKinkDaughters(kFALSE),
 164 dCutMaxChi2TPCConstrainedGlobal(36.),
 165 // binning for THnSparse
 166 fMultNbins(0),
 167 fPtNbins(0),
 168 fPtCorrNbins(0),
 169 fEtaNbins(0),
 170 fZvNbins(0),
 171 fCentralityNbins(0),
 172 fBinsMult(0),
 173 fBinsPt(0),
 174 fBinsPtCorr(0),
 175 fBinsEta(0),
 176 fBinsZv(0),
 177 fBinsCentrality(0)
 178 {
 179   fMultNbins = 0;
 180   fPtNbins = 0;
 181   fPtCorrNbins = 0;
 182   fEtaNbins = 0;
 183   fZvNbins = 0;
 184   fCentralityNbins = 0;
 185   fBinsMult = 0;
 186   fBinsPt = 0;
 187   fBinsPtCorr = 0;
 188   fBinsEta = 0;
 189   fBinsEta = 0;
 190   fBinsZv = 0;
 191   fBinsCentrality = 0;
 192
 193   DefineOutput(1, TList::Class());
 194 }
 195
 196 // destructor
 197 AlidNdPtAnalysisPbPbAOD::~AlidNdPtAnalysisPbPbAOD()
 198 {
 199   if(hnZvPtEtaCent) delete hnZvPtEtaCent; hnZvPtEtaCent = 0;
 200   if(hPt) delete hPt; hPt = 0;
 201   if(hnMCRecPrimZvPtEtaCent) delete hnMCRecPrimZvPtEtaCent; hnMCRecPrimZvPtEtaCent = 0;
 202   if(hnMCGenZvPtEtaCent) delete hnMCGenZvPtEtaCent; hnMCGenZvPtEtaCent = 0;
 203   if(hnMCRecSecZvPtEtaCent) delete hnMCRecSecZvPtEtaCent; hnMCRecSecZvPtEtaCent = 0;
 204   if(hMCPt) delete hMCPt; hMCPt = 0;
 205   if(hEventStatistics) delete hEventStatistics; hEventStatistics = 0;
 206   if(hEventStatisticsCentrality) delete hEventStatisticsCentrality; hEventStatisticsCentrality = 0;
 207   if(hAllEventStatisticsCentrality) delete hAllEventStatisticsCentrality; hAllEventStatisticsCentrality = 0;
 208   if(hnZvMultCent) delete hnZvMultCent; hnZvMultCent = 0;
 209   if(hTriggerStatistics) delete hTriggerStatistics; hTriggerStatistics = 0;
 210   if(hMCTrackPdgCode) delete hMCTrackPdgCode; hMCTrackPdgCode = 0;
 211   if(hMCTrackStatusCode) delete hMCTrackStatusCode; hMCTrackStatusCode = 0;
 212   if(hCharge) delete hCharge; hCharge = 0;
 213   if(hMCCharge) delete hMCCharge; hMCCharge = 0;
 214   if(hMCPdgPt) delete hMCPdgPt; hMCPdgPt = 0;
 215   if(hMCHijingPrim) delete hMCHijingPrim; hMCHijingPrim = 0;
 216   if(hAccNclsTPC) delete hAccNclsTPC; hAccNclsTPC = 0;
 217   if(hAccCrossedRowsTPC) delete hAccCrossedRowsTPC; hAccCrossedRowsTPC = 0;
 218   if(hDCAPtAll) delete hDCAPtAll; hDCAPtAll = 0;
 219   if(hDCAPtAccepted) delete hDCAPtAccepted; hDCAPtAccepted = 0;
 220   if(hMCDCAPtSecondary) delete hMCDCAPtSecondary; hMCDCAPtSecondary = 0;
 221   if(hMCDCAPtPrimary) delete hMCDCAPtPrimary; hMCDCAPtPrimary = 0;
 222 }
 223
 224 void AlidNdPtAnalysisPbPbAOD::UserCreateOutputObjects()
 225 {
 226   // create all output histograms here
 227   OpenFile(1, "RECREATE");
 228
 229   fOutputList = new TList();
 230   fOutputList->SetOwner();
 231
 232   //define default binning
 233   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 };
 234   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};
 235   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};
 236   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};
 237   Double_t binsZvDefault[13] = {-30.,-25.,-20.,-15.,-10.,-5.,0.,5.,10.,15.,20.,25.,30.};
 238   Double_t binsCentralityDefault[12] = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.};
 239
 240   // if no binning is set, use the default
 241   if (!fBinsMult)       { SetBinsMult(48,binsMultDefault); }
 242   if (!fBinsPt)         { SetBinsPt(82,binsPtDefault); }
 243   if (!fBinsPtCorr)     { SetBinsPtCorr(37,binsPtCorrDefault); }
 244   if (!fBinsEta)        { SetBinsEta(31,binsEtaDefault); }
 245   if (!fBinsZv)         { SetBinsZv(13,binsZvDefault); }
 246   if (!fBinsCentrality) { SetBinsCentrality(12,binsCentralityDefault); }
 247
 248   Int_t binsZvPtEtaCent[4]={fZvNbins-1,fPtNbins-1,fEtaNbins-1,fCentralityNbins-1};
 249   Int_t binsZvMultCent[3]={fZvNbins-1,fMultNbins-1,fCentralityNbins-1};
 250
 251   // define Histograms
 252   hnZvPtEtaCent = new THnSparseF("hnZvPtEtaCent","Zv:Pt:Eta:Centrality",4,binsZvPtEtaCent);
 253   hnZvPtEtaCent->SetBinEdges(0,fBinsZv);
 254   hnZvPtEtaCent->SetBinEdges(1,fBinsPt);
 255   hnZvPtEtaCent->SetBinEdges(2,fBinsEta);
 256   hnZvPtEtaCent->SetBinEdges(3,fBinsCentrality);
 257   hnZvPtEtaCent->GetAxis(0)->SetTitle("Zv (cm)");
 258   hnZvPtEtaCent->GetAxis(1)->SetTitle("Pt (GeV/c)");
 259   hnZvPtEtaCent->GetAxis(2)->SetTitle("Eta");
 260   hnZvPtEtaCent->GetAxis(3)->SetTitle("Centrality");
 261   hnZvPtEtaCent->Sumw2();
 262
 263   hnMCRecPrimZvPtEtaCent = new THnSparseF("hnMCRecPrimZvPtEtaCent","mcZv:mcPt:mcEta:Centrality",4,binsZvPtEtaCent);
 264   hnMCRecPrimZvPtEtaCent->SetBinEdges(0,fBinsZv);
 265   hnMCRecPrimZvPtEtaCent->SetBinEdges(1,fBinsPt);
 266   hnMCRecPrimZvPtEtaCent->SetBinEdges(2,fBinsEta);
 267   hnMCRecPrimZvPtEtaCent->SetBinEdges(3,fBinsCentrality);
 268   hnMCRecPrimZvPtEtaCent->GetAxis(0)->SetTitle("MC Zv (cm)");
 269   hnMCRecPrimZvPtEtaCent->GetAxis(1)->SetTitle("MC Pt (GeV/c)");
 270   hnMCRecPrimZvPtEtaCent->GetAxis(2)->SetTitle("MC Eta");
 271   hnMCRecPrimZvPtEtaCent->GetAxis(3)->SetTitle("Centrality");
 272   hnMCRecPrimZvPtEtaCent->Sumw2();
 273
 274   hnMCGenZvPtEtaCent = new THnSparseF("hnMCGenZvPtEtaCent","mcZv:mcPt:mcEta:Centrality",4,binsZvPtEtaCent);
 275   hnMCGenZvPtEtaCent->SetBinEdges(0,fBinsZv);
 276   hnMCGenZvPtEtaCent->SetBinEdges(1,fBinsPt);
 277   hnMCGenZvPtEtaCent->SetBinEdges(2,fBinsEta);
 278   hnMCGenZvPtEtaCent->SetBinEdges(3,fBinsCentrality);
 279   hnMCGenZvPtEtaCent->GetAxis(0)->SetTitle("MC Zv (cm)");
 280   hnMCGenZvPtEtaCent->GetAxis(1)->SetTitle("MC Pt (GeV/c)");
 281   hnMCGenZvPtEtaCent->GetAxis(2)->SetTitle("MC Eta");
 282   hnMCGenZvPtEtaCent->GetAxis(3)->SetTitle("Centrality");
 283   hnMCGenZvPtEtaCent->Sumw2();
 284
 285   hnMCRecSecZvPtEtaCent = new THnSparseF("hnMCRecSecZvPtEtaCent","mcZv:mcPt:mcEta:Centrality",4,binsZvPtEtaCent);
 286   hnMCRecSecZvPtEtaCent->SetBinEdges(0,fBinsZv);
 287   hnMCRecSecZvPtEtaCent->SetBinEdges(1,fBinsPt);
 288   hnMCRecSecZvPtEtaCent->SetBinEdges(2,fBinsEta);
 289   hnMCRecSecZvPtEtaCent->SetBinEdges(3,fBinsCentrality);
 290   hnMCRecSecZvPtEtaCent->GetAxis(0)->SetTitle("MC Sec Zv (cm)");
 291   hnMCRecSecZvPtEtaCent->GetAxis(1)->SetTitle("MC Sec Pt (GeV/c)");
 292   hnMCRecSecZvPtEtaCent->GetAxis(2)->SetTitle("MC Sec Eta");
 293   hnMCRecSecZvPtEtaCent->GetAxis(3)->SetTitle("Centrality");
 294   hnMCRecSecZvPtEtaCent->Sumw2();
 295
 296   hPt = new TH1F("hPt","hPt",2000,0,200);
 297   hPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 298   hPt->GetYaxis()->SetTitle("dN/dp_{T}");
 299   hPt->Sumw2();
 300
 301   hMCPt = new TH1F("hMCPt","hMCPt",2000,0,200);
 302   hMCPt->GetXaxis()->SetTitle("MC p_{T} (GeV/c)");
 303   hMCPt->GetYaxis()->SetTitle("dN/dp_{T}");
 304   hMCPt->Sumw2();
 305
 306   hEventStatistics = new TH1F("hEventStatistics","hEventStatistics",10,0,10);
 307   hEventStatistics->GetYaxis()->SetTitle("number of events");
 308   hEventStatistics->SetBit(TH1::kCanRebin);
 309
 310   hEventStatisticsCentrality = new TH1F("hEventStatisticsCentrality","hEventStatisticsCentrality",fCentralityNbins-1, fBinsCentrality);
 311   hEventStatisticsCentrality->GetYaxis()->SetTitle("number of events");
 312
 313   hAllEventStatisticsCentrality = new TH1F("hAllEventStatisticsCentrality","hAllEventStatisticsCentrality",fCentralityNbins-1, fBinsCentrality);
 314   hAllEventStatisticsCentrality->GetYaxis()->SetTitle("number of events");
 315
 316   hnZvMultCent = new THnSparseF("hnZvMultCent","Zv:mult:Centrality",3,binsZvMultCent);
 317   hnZvMultCent->SetBinEdges(0,fBinsZv);
 318   hnZvMultCent->SetBinEdges(1,fBinsMult);
 319   hnZvMultCent->SetBinEdges(2,fBinsCentrality);
 320   hnZvMultCent->GetAxis(0)->SetTitle("Zv (cm)");
 321   hnZvMultCent->GetAxis(1)->SetTitle("N_{acc}");
 322   hnZvMultCent->GetAxis(2)->SetTitle("Centrality");
 323   hnZvMultCent->Sumw2();
 324
 325   hTriggerStatistics = new TH1F("hTriggerStatistics","hTriggerStatistics",10,0,10);
 326   hTriggerStatistics->GetYaxis()->SetTitle("number of events");
 327
 328   hMCTrackPdgCode = new TH1F("hMCTrackPdgCode","hMCTrackPdgCode",100,0,10);
 329   hMCTrackPdgCode->GetYaxis()->SetTitle("number of tracks");
 330   hMCTrackPdgCode->SetBit(TH1::kCanRebin);
 331
 332   hMCTrackStatusCode = new TH1F("hMCTrackStatusCode","hMCTrackStatusCode",100,0,10);
 333   hMCTrackStatusCode->GetYaxis()->SetTitle("number of tracks");
 334   hMCTrackStatusCode->SetBit(TH1::kCanRebin);
 335
 336   hCharge = new TH1F("hCharge","hCharge",30, -5, 5);
 337   hCharge->GetXaxis()->SetTitle("Charge");
 338   hCharge->GetYaxis()->SetTitle("number of tracks");
 339
 340   hMCCharge = new TH1F("hMCCharge","hMCCharge",30, -5, 5);
 341   hMCCharge->GetXaxis()->SetTitle("MC Charge");
 342   hMCCharge->GetYaxis()->SetTitle("number of tracks");
 343
 344   hMCPdgPt = new TH2F("hMCPdgPt","hMCPdgPt",fPtNbins-1, fBinsPt, 100,0,100);
 345   hMCPdgPt->GetYaxis()->SetTitle("particle");
 346   hMCPdgPt->GetXaxis()->SetTitle("Pt (GeV/c)");
 347
 348   hMCHijingPrim = new TH1F("hMCHijingPrim","hMCHijingPrim",2,0,2);
 349   hMCPdgPt->GetYaxis()->SetTitle("number of particles");
 350
 351   hAccNclsTPC = new TH1F("hAccNclsTPC","hAccNclsTPC",160,0,159);
 352   hAccNclsTPC->GetXaxis()->SetTitle("number of clusters per track after cut");
 353
 354   hAccCrossedRowsTPC = new TH1F("hAccCrossedRowsTPC","hAccCrossedRowsTPC",160,0,159);
 355   hAccCrossedRowsTPC->GetXaxis()->SetTitle("number of crossed rows per track after cut");
 356
 357   Int_t binsDCAxyDCAzPt[3] = { 200,200, fPtNbins-1};
 358   Double_t minDCAxyDCAzPt[3] = { -10, -10, 0};
 359   Double_t maxDCAxyDCAzPt[3] = { 10., 10., 100};
 360
 361   hDCAPtAll = new THnSparseF("hDCAPtAll","hDCAPtAll;p_{T} (GeV/c);DCAz",3, binsDCAxyDCAzPt, minDCAxyDCAzPt, maxDCAxyDCAzPt);
 362   hDCAPtAccepted = new THnSparseF("hDCAPtAccepted","hDCAPtAccepted;p_{T} (GeV/c);DCAz",3, binsDCAxyDCAzPt, minDCAxyDCAzPt, maxDCAxyDCAzPt);
 363   hMCDCAPtSecondary = new THnSparseF("hMCDCAPtSecondary","hMCDCAPtSecondary;p_{T} (GeV/c);DCAz",3, binsDCAxyDCAzPt, minDCAxyDCAzPt, maxDCAxyDCAzPt);
 364   hMCDCAPtPrimary = new THnSparseF("hMCDCAPtPrimary","hMCDCAPtPrimary;p_{T} (GeV/c);DCAz",3, binsDCAxyDCAzPt, minDCAxyDCAzPt, maxDCAxyDCAzPt);
 365
 366   hDCAPtAll->SetBinEdges(2, fBinsPt);
 367   hDCAPtAccepted->SetBinEdges(2, fBinsPt);
 368   hMCDCAPtSecondary->SetBinEdges(2, fBinsPt);
 369   hMCDCAPtPrimary->SetBinEdges(2, fBinsPt);
 370
 371   // Add Histos, Profiles etc to List
 372   fOutputList->Add(hnZvPtEtaCent);
 373   fOutputList->Add(hPt);
 374   fOutputList->Add(hnMCRecPrimZvPtEtaCent);
 375   fOutputList->Add(hnMCGenZvPtEtaCent);
 376   fOutputList->Add(hnMCRecSecZvPtEtaCent);
 377   fOutputList->Add(hMCPt);
 378   fOutputList->Add(hEventStatistics);
 379   fOutputList->Add(hEventStatisticsCentrality);
 380   fOutputList->Add(hAllEventStatisticsCentrality);
 381   fOutputList->Add(hnZvMultCent);
 382   fOutputList->Add(hTriggerStatistics);
 383   fOutputList->Add(hMCTrackPdgCode);
 384   fOutputList->Add(hMCTrackStatusCode);
 385   fOutputList->Add(hCharge);
 386   fOutputList->Add(hMCCharge);
 387   fOutputList->Add(hMCPdgPt);
 388   fOutputList->Add(hMCHijingPrim);
 389   fOutputList->Add(hAccNclsTPC);
 390   fOutputList->Add(hAccCrossedRowsTPC);
 391   fOutputList->Add(hDCAPtAll);
 392   fOutputList->Add(hDCAPtAccepted);
 393   fOutputList->Add(hMCDCAPtSecondary);
 394   fOutputList->Add(hMCDCAPtPrimary);
 395
 396
 397   PostData(1, fOutputList);
 398 }
 399
 400 void AlidNdPtAnalysisPbPbAOD::UserExec(Option_t *option)
 401 {
 402
 403   // Main Loop
 404   // called for each event
 405   hEventStatistics->Fill("all events",1);
 406
 407   // set ZERO pointers:
 408   AliInputEventHandler *inputHandler = NULL;
 409   AliAODTrack *track = NULL;
 410   AliAODMCParticle *mcPart = NULL;
 411   AliAODMCHeader *mcHdr = NULL;
 412   AliGenHijingEventHeader *genHijingHeader = NULL;
 413   AliGenPythiaEventHeader *genPythiaHeader = NULL;
 414
 415   Bool_t bIsEventSelectedMB = kFALSE;
 416   Bool_t bIsEventSelectedSemi = kFALSE;
 417   Bool_t bIsEventSelectedCentral = kFALSE;
 418   Bool_t bIsEventSelected = kFALSE;
 419   Bool_t bIsPrimary = kFALSE;
 420   Bool_t bIsHijingParticle = kFALSE;
 421   Bool_t bIsPythiaParticle = kFALSE;
 422   Bool_t bEventHasATrack = kFALSE;
 423   Bool_t bEventHasATrackInRange = kFALSE;
 424   Int_t nTriggerFired = 0;
 425
 426
 427   Double_t dMCTrackZvPtEtaCent[4] = {0};
 428   Double_t dTrackZvPtEtaCent[4] = {0};
 429
 430   Double_t dMCEventZv = -100;
 431   Double_t dEventZv = -100;
 432   Int_t iAcceptedMultiplicity = 0;
 433
 434   bIsMonteCarlo = kFALSE;
 435
 436   AliAODEvent *eventAOD = 0x0;
 437   eventAOD = dynamic_cast<AliAODEvent*>( InputEvent() );
 438   if (!eventAOD) {
 439     AliWarning("ERROR: eventAOD not available \n");
 440     return;
 441   }
 442
 443   // check, which trigger has been fired
 444   inputHandler = (AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 445   bIsEventSelectedMB = ( inputHandler->IsEventSelected() & AliVEvent::kMB);
 446   bIsEventSelectedSemi = ( inputHandler->IsEventSelected() & AliVEvent::kSemiCentral);
 447   bIsEventSelectedCentral = ( inputHandler->IsEventSelected() & AliVEvent::kCentral);
 448
 449   if(bIsEventSelectedMB || bIsEventSelectedSemi || bIsEventSelectedCentral) hTriggerStatistics->Fill("all triggered events",1);
 450   if(bIsEventSelectedMB) { hTriggerStatistics->Fill("MB trigger",1); nTriggerFired++; }
 451   if(bIsEventSelectedSemi) { hTriggerStatistics->Fill("SemiCentral trigger",1); nTriggerFired++; }
 452   if(bIsEventSelectedCentral) { hTriggerStatistics->Fill("Central trigger",1); nTriggerFired++; }
 453   if(nTriggerFired == 0) { hTriggerStatistics->Fill("No trigger",1); }
 454
 455   bIsEventSelected = ( inputHandler->IsEventSelected() & GetCollisionCandidates() );
 456
 457   // only take tracks of events, which are triggered
 458   if(nTriggerFired == 0) { return; }
 459
 460   //   if( !bIsEventSelected || nTriggerFired>1 ) return;
 461
 462   //   hEventStatistics->Fill("events with only coll. cand.", 1);
 463
 464
 465
 466   // check if there is a stack, if yes, then do MC loop
 467   TList *list = eventAOD->GetList();
 468   TClonesArray *stack = 0x0;
 469   stack = (TClonesArray*)list->FindObject(AliAODMCParticle::StdBranchName());
 470
 471   if( stack )
 472   {
 473     bIsMonteCarlo = kTRUE;
 474
 475     mcHdr = (AliAODMCHeader*)list->FindObject(AliAODMCHeader::StdBranchName());
 476
 477     genHijingHeader = GetHijingEventHeader(mcHdr);
 478     //     genPythiaHeader = GetPythiaEventHeader(mcHdr);
 479
 480     if(!genHijingHeader) { return; }
 481
 482     //     if(!genPythiaHeader)  { return; }
 483
 484     dMCEventZv = mcHdr->GetVtxZ();
 485     dMCTrackZvPtEtaCent[0] = dMCEventZv;
 486     hEventStatistics->Fill("MC all events",1);
 487   }
 488
 489   AliCentrality* aCentrality = eventAOD->GetCentrality();
 490   Double_t dCentrality = aCentrality->GetCentralityPercentile("V0M");
 491
 492   if( dCentrality < 0 ) return;
 493   hEventStatistics->Fill("after centrality selection",1);
 494
 495
 496
 497   // start with MC truth analysis
 498   if(bIsMonteCarlo)
 499   {
 500
 501     if( dMCEventZv > dCutMaxZVertex )  { return; }
 502
 503     dMCTrackZvPtEtaCent[0] = dMCEventZv;
 504
 505     hEventStatistics->Fill("MC afterZv cut",1);
 506
 507     for(Int_t iMCtrack = 0; iMCtrack < stack->GetEntriesFast(); iMCtrack++)
 508     {
 509       mcPart =(AliAODMCParticle*)stack->At(iMCtrack);
 510
 511       // check for charge
 512       if( !(IsMCTrackAccepted(mcPart)) ) continue;
 513
 514       if(!IsHijingParticle(mcPart, genHijingHeader)) { continue; }
 515
 516       if(mcPart->IsPhysicalPrimary() )
 517       {
 518         hMCHijingPrim->Fill("IsPhysicalPrimary",1);
 519       }
 520       else
 521       {
 522         hMCHijingPrim->Fill("NOT a primary",1);
 523         continue;
 524       }
 525
 526
 527       //
 528       // ======================== fill histograms ========================
 529       dMCTrackZvPtEtaCent[1] = mcPart->Pt();
 530       dMCTrackZvPtEtaCent[2] = mcPart->Eta();
 531       dMCTrackZvPtEtaCent[3] = dCentrality;
 532
 533       bEventHasATrack = kTRUE;
 534
 535       hnMCGenZvPtEtaCent->Fill(dMCTrackZvPtEtaCent);
 536
 537       if( (dMCTrackZvPtEtaCent[1] > dCutPtMin) &&
 538         (dMCTrackZvPtEtaCent[1] < dCutPtMax) &&
 539         (dMCTrackZvPtEtaCent[2] > dCutEtaMin) &&
 540         (dMCTrackZvPtEtaCent[2] < dCutEtaMax) )
 541       {
 542         hMCPt->Fill(mcPart->Pt());
 543         hMCCharge->Fill(mcPart->Charge()/3.);
 544         bEventHasATrackInRange = kTRUE;
 545       }
 546
 547     }
 548   } // isMonteCarlo
 549   if(bEventHasATrack) { hEventStatistics->Fill("MC events with tracks",1); }
 550   if(bEventHasATrackInRange) { hEventStatistics->Fill("MC events with tracks in range",1); }
 551   bEventHasATrack = kFALSE;
 552   bEventHasATrackInRange = kFALSE;
 553
 554
 555
 556   // Loop over recontructed tracks
 557
 558   dEventZv = eventAOD->GetPrimaryVertex()->GetZ();
 559   if( TMath::Abs(dEventZv) > dCutMaxZVertex ) return;
 560
 561   hAllEventStatisticsCentrality->Fill(dCentrality/*, nTriggerFired*/);
 562
 563   hEventStatistics->Fill("after Zv cut",1);
 564
 565   dTrackZvPtEtaCent[0] = dEventZv;
 566
 567   for(Int_t itrack = 0; itrack < eventAOD->GetNumberOfTracks(); itrack++)
 568   {
 569     track = eventAOD->GetTrack(itrack);
 570     if(!track) continue;
 571
 572     mcPart = NULL;
 573     dMCTrackZvPtEtaCent[1] = 0;
 574     dMCTrackZvPtEtaCent[2] = 0;
 575     dMCTrackZvPtEtaCent[3] = 0;
 576
 577     bIsPrimary = kFALSE;
 578
 579     Double_t dDCAxyDCAzPt[3] = { GetDCAxy(track, eventAOD), GetDCAz(track, eventAOD), track->Pt() };
 580
 581     hDCAPtAll->Fill(dDCAxyDCAzPt);
 582
 583     if( !(IsTrackAccepted(track)) ) continue;
 584
 585     dTrackZvPtEtaCent[1] = track->Pt();
 586     dTrackZvPtEtaCent[2] = track->Eta();
 587     dTrackZvPtEtaCent[3] = dCentrality;
 588
 589     if( bIsMonteCarlo )
 590     {
 591       mcPart = (AliAODMCParticle*)stack->At(TMath::Abs(track->GetLabel()));
 592       if( !mcPart ) { continue; }
 593
 594       // check for charge
 595       if( !(IsMCTrackAccepted(mcPart)) ) {  continue; }
 596
 597       bIsHijingParticle = IsHijingParticle(mcPart, genHijingHeader);
 598       //       bIsPythiaParticle = IsPythiaParticle(mcPart, genPythiaHeader);
 599
 600       //       if(!bIsHijingParticle) continue; // only take real tracks, not injected ones
 601
 602       bIsPrimary = mcPart->IsPhysicalPrimary();
 603
 604       dMCTrackZvPtEtaCent[1] = mcPart->Pt();
 605       dMCTrackZvPtEtaCent[2] = mcPart->Eta();
 606       dMCTrackZvPtEtaCent[3] = dCentrality;
 607
 608       if(bIsPrimary && bIsHijingParticle)
 609       {
 610         hnMCRecPrimZvPtEtaCent->Fill(dMCTrackZvPtEtaCent);
 611         hMCDCAPtPrimary->Fill(dDCAxyDCAzPt);
 612       }
 613
 614       if(!bIsPrimary /*&& !bIsHijingParticle*/)
 615       {
 616         Int_t indexMoth = mcPart->GetMother();
 617         if(indexMoth >= 0)
 618         {
 619           AliAODMCParticle* moth = (AliAODMCParticle*)stack->At(indexMoth);
 620           Bool_t bMotherIsHijingParticle = IsHijingParticle(moth, genHijingHeader);
 621
 622           if(bMotherIsHijingParticle) // only store secondaries, which come from a not embedded signal!
 623           {
 624             hMCTrackStatusCode->Fill(Form("%d",mcPart->GetStatus()), 1);
 625             if(TMath::Abs(mcPart->Eta()) < 0.8) { hMCPdgPt->Fill(mcPart->Pt(), Form("%s",GetParticleName(mcPart->GetPdgCode())), 1); }
 626
 627             hnMCRecSecZvPtEtaCent->Fill(dMCTrackZvPtEtaCent);
 628             hMCDCAPtSecondary->Fill(dDCAxyDCAzPt);
 629             hMCTrackPdgCode->Fill(Form("%s_H%i_H%i",GetParticleName(moth->GetPdgCode()),bMotherIsHijingParticle, bIsHijingParticle), 1);
 630             //    delete moth;
 631           }
 632         }
 633       }
 634     } // end isMonteCarlo
 635
 636     // ======================== fill histograms ========================
 637
 638     //     if(bIsMonteCarlo && !bIsHijingParticle)
 639     //     {
 640       //       continue;  //only store reco tracks, which do not come from embedded signal
 641       //     }
 642
 643
 644
 645       bEventHasATrack = kTRUE;
 646
 647       hnZvPtEtaCent->Fill(dTrackZvPtEtaCent);
 648       hDCAPtAccepted->Fill(dDCAxyDCAzPt);
 649
 650       if( (dTrackZvPtEtaCent[1] > dCutPtMin) &&
 651         (dTrackZvPtEtaCent[1] < dCutPtMax) &&
 652         (dTrackZvPtEtaCent[2] > dCutEtaMin) &&
 653         (dTrackZvPtEtaCent[2] < dCutEtaMax) )
 654       {
 655         iAcceptedMultiplicity++;
 656         bEventHasATrackInRange = kTRUE;
 657         hPt->Fill(track->Pt());
 658         hCharge->Fill(track->Charge());
 659       }
 660   } // end track loop
 661
 662   if(bEventHasATrack) { hEventStatistics->Fill("events with tracks",1); bEventHasATrack = kFALSE; }
 663
 664   if(bEventHasATrackInRange)
 665   {
 666     hEventStatistics->Fill("events with tracks in range",1);
 667     hEventStatisticsCentrality->Fill(dCentrality);
 668     bEventHasATrackInRange = kFALSE;
 669   }
 670
 671   Double_t dEventZvMultCent[3] = {dEventZv, iAcceptedMultiplicity, dCentrality};
 672   hnZvMultCent->Fill(dEventZvMultCent);
 673
 674
 675
 676   PostData(1, fOutputList);
 677
 678 }
 679
 680 Bool_t AlidNdPtAnalysisPbPbAOD::IsTrackAccepted(AliAODTrack *tr)
 681 {
 682   if(!tr) return kFALSE;
 683
 684   if(tr->Charge()==0) { return kFALSE; }
 685
 686   if(!(tr->TestFilterBit(AliAODTrack::kTrkGlobal)) ) { return kFALSE; }
 687
 688   Double_t dNClustersTPC = tr->GetTPCNcls();
 689   Double_t dCrossedRowsTPC = tr->GetTPCClusterInfo(2,1);
 690
 691   if(dCrossedRowsTPC < GetCutMinNCrossedRowsTPC()) { return kFALSE; }
 692
 693   hAccNclsTPC->Fill(dNClustersTPC);
 694   hAccCrossedRowsTPC->Fill(dCrossedRowsTPC);
 695   //   Double_t dFindableClustersTPC = tr->GetTPCNclsF();
 696   //   Double_t dChi2PerClusterTPC = (dNClustersTPC>0)?tr->Chi2perNDF()*(dNClustersTPC-5)/dNClustersTPC:-1.; // see AliDielectronVarManager.h
 697   //
 698   //   Bool_t bIsFromKink = kFALSE;
 699   //   if(tr->GetProdVertex()->GetType() == AliAODVertex::kKink) bIsFromKink = kTRUE;
 700   //
 701   //   // from AliAnalysisTaskPIDqa.cxx
 702   //   ULong_t uStatus = tr->GetStatus();
 703   //   Bool_t bHasRefitTPC = kFALSE;
 704   //   Bool_t bHasRefitITS = kFALSE;
 705   //
 706   //   if ((uStatus & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) bHasRefitTPC = kTRUE;
 707   //   if ((uStatus & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) bHasRefitITS = kTRUE;
 708   //
 709   //   // from AliDielectronVarManager.h
 710   //   Bool_t bHasHitInSPD = kFALSE;
 711   //   for (Int_t iC=0; iC<2; iC++)
 712   //   {
 713     //     if (((tr->GetITSClusterMap()) & (1<<(iC))) > 0) {  bHasHitInSPD = kTRUE;  }
 714     //   }
 715     //
 716     //   Double_t dNClustersITS = tr->GetITSNcls();
 717
 718     // cuts to be done:
 719     // TPC
 720     //   esdTrackCuts->SetMinNCrossedRowsTPC(70);
 721     //   esdTrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
 722     //
 723     //   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
 724     //   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 725     //   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 726     // ITS
 727     //   esdTrackCuts->SetRequireITSRefit(kTRUE);
 728     //   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,   AliESDtrackCuts::kAny);
 729     //
 730     //   esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0105+0.0350/pt^1.1");
 731     //   esdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
 732     //
 733     //   esdTrackCuts->SetMaxDCAToVertexZ(2);
 734     //   esdTrackCuts->SetDCAToVertex2D(kFALSE);
 735     //   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 736     //
 737     //   esdTrackCuts->SetMaxChi2PerClusterITS(36);
 738
 739
 740     return kTRUE;
 741 }
 742
 743 Double_t AlidNdPtAnalysisPbPbAOD::GetDCAz(AliAODTrack *track, AliAODEvent *event)
 744 {
 745  return GetDCA(track, event, kTRUE);
 746 }
 747
 748 Double_t AlidNdPtAnalysisPbPbAOD::GetDCAxy(AliAODTrack *track, AliAODEvent *event)
 749 {
 750  return GetDCA(track, event, kFALSE);
 751 }
 752
 753 Double_t AlidNdPtAnalysisPbPbAOD::GetDCA(AliAODTrack *tr, AliAODEvent *evt, Bool_t bDCAz)
 754 {
 755   if(!tr) return -999.;
 756
 757   if(tr->TestBit(AliAODTrack::kIsDCA))
 758   {
 759     if(bDCAz) return tr->ZAtDCA();
 760     else return sqrt(tr->XAtDCA()*tr->XAtDCA() + tr->YAtDCA()*tr->YAtDCA());
 761   }
 762
 763   Bool_t ok=kFALSE;
 764   Double_t d0z0[2];
 765   if(evt)
 766   {
 767     Double_t covd0z0[3];
 768
 769     AliExternalTrackParam etp;
 770     etp.CopyFromVTrack(tr);
 771
 772     Float_t xstart = etp.GetX();
 773     if(xstart>3.)
 774     {
 775       d0z0[0]=-999.;
 776       d0z0[1]=-999.;
 777       //printf("This method can be used only for propagation inside the beam pipe \n");
 778       if(bDCAz) return d0z0[1];
 779       else return d0z0[0];
 780     }
 781
 782
 783     AliAODVertex *vtx =(AliAODVertex*)(evt->GetPrimaryVertex());
 784     Double_t fBzkG = evt->GetMagneticField(); // z componenent of field in kG
 785     ok = etp.PropagateToDCA(vtx,fBzkG,kVeryBig,d0z0,covd0z0);
 786   }
 787   if(!ok){
 788     d0z0[0]=-999.;
 789     d0z0[1]=-999.;
 790     if(bDCAz) return d0z0[1];
 791     else return d0z0[0];
 792   }
 793   if(bDCAz) return d0z0[1];
 794   else return d0z0[0];
 795 }
 796
 797 Bool_t AlidNdPtAnalysisPbPbAOD::IsMCTrackAccepted(AliAODMCParticle *part)
 798 {
 799   if(!part) return kFALSE;
 800
 801   Double_t charge = part->Charge()/3.;
 802   if (TMath::Abs(charge) < 0.001) return kFALSE;
 803
 804   return kTRUE;
 805 }
 806
 807 const char * AlidNdPtAnalysisPbPbAOD::GetParticleName(Int_t pdg)
 808 {
 809   TParticlePDG * p1 = TDatabasePDG::Instance()->GetParticle(pdg);
 810   if(p1) return p1->GetName();
 811   return Form("%d", pdg);
 812 }
 813
 814 AliGenHijingEventHeader* AlidNdPtAnalysisPbPbAOD::GetHijingEventHeader(AliAODMCHeader *header)
 815 {
 816   //
 817   // inspired by PWGJE/AliPWG4HighPtSpectra.cxx
 818   //
 819
 820   if(!header) return 0x0;
 821   AliGenHijingEventHeader* hijingGenHeader = NULL;
 822
 823   TList* headerList = header->GetCocktailHeaders();
 824
 825   for(Int_t i = 0; i < headerList->GetEntries(); i++)
 826   {
 827     hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(headerList->At(i));
 828     if(hijingGenHeader) break;
 829   }
 830
 831   if(!hijingGenHeader) return 0x0;
 832
 833   return hijingGenHeader;
 834 }
 835
 836 AliGenPythiaEventHeader* AlidNdPtAnalysisPbPbAOD::GetPythiaEventHeader(AliAODMCHeader *header)
 837 {
 838   //
 839   // inspired by PWGJE/AliPWG4HighPtSpectra.cxx
 840   //
 841
 842   if(!header) return 0x0;
 843   AliGenPythiaEventHeader* PythiaGenHeader = NULL;
 844
 845   TList* headerList = header->GetCocktailHeaders();
 846
 847   for(Int_t i = 0; i < headerList->GetEntries(); i++)
 848   {
 849     PythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
 850     if(PythiaGenHeader) break;
 851   }
 852
 853   if(!PythiaGenHeader) return 0x0;
 854
 855   return PythiaGenHeader;
 856 }
 857
 858 //________________________________________________________________________
 859 Bool_t AlidNdPtAnalysisPbPbAOD::IsHijingParticle(const AliAODMCParticle *part, AliGenHijingEventHeader* hijingGenHeader){
 860
 861   // Check whether a particle is from Hijing or some injected
 862
 863   if(part->Label() > (hijingGenHeader->NProduced()-1)) return kFALSE;
 864   return kTRUE;
 865 }
 866
 867 //________________________________________________________________________
 868 Bool_t AlidNdPtAnalysisPbPbAOD::IsPythiaParticle(const AliAODMCParticle *part, AliGenPythiaEventHeader* pythiaGenHeader){
 869
 870   // Check whether a particle is from Pythia or some injected
 871
 872   if(part->Label() > (pythiaGenHeader->NProduced()-1)) return kFALSE;
 873   return kTRUE;
 874 }
 875
 876 // Int_t AlidNdPtAnalysisPbPbAOD::IsMCSecondary(AliAODMCParticle *part, TClonesArray *arrayMC)
 877 // {
 878   //   //
 879   //   // adapted from AliAnalysisTaskSpectraAOD.cxx
 880   //   //
 881   //   // returns
 882   //   // -1: no particle
 883   //   // 0: is primary
 884   //   // 1: is secondary from weak
 885   //   // 2: is secondary from material
 886   //
 887   //   // usage for studies, currrently not implemented
 888   //
 889   //   if(!part) return -1;
 890   //
 891   //   if( part->IsPhysicalPrimary() ) return 0;
 892   //
 893   //   Bool_t isSecondaryMaterial = kFALSE;
 894   //   Bool_t isSecondaryWeak     = kFALSE;
 895   //   Int_t mfl = -999;
 896   //   Int_t codemoth = -999;
 897   //   Int_t indexMoth = part->GetMother(); // FIXME ignore fakes? TO BE CHECKED, on ESD is GetFirstMother()
 898   //   if(indexMoth >= 0)
 899   //   {
 900     //     AliAODMCParticle* moth = (AliAODMCParticle*) arrayMC->At(indexMoth);
 901     //     codemoth = TMath::Abs(moth->GetPdgCode());
 902     //     mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
 903     //   }
 904     //   // add if(partMC->GetStatus() & kPDecay)? FIXME
 905     //   if(mfl==3) isSecondaryWeak     = kTRUE;
 906     //   else       isSecondaryMaterial = kTRUE;
 907     //
 908     //   if(isSecondaryWeak) return 1;
 909     //   if(isSecondaryMaterial) return 2;
 910     //
 911     // //   if( isSecondaryMaterial || isSecondaryWeak ) return kTRUE;
 912     //
 913     //   // return kFALSE; this line will not be reached, as either isSecondaryMaterial or isSecondaryWeak is true!
 914     //   // removed due to coverity
 915     // }
 916
 917
 918
 919     void AlidNdPtAnalysisPbPbAOD::Terminate(Option_t *)
 920     {
 921
 922     }
 923
 924     Double_t* AlidNdPtAnalysisPbPbAOD::GetArrayClone(Int_t n, Double_t* source)
 925     {
 926       if (!source || n==0) return 0;
 927       Double_t* dest = new Double_t[n];
 928       for (Int_t i=0; i<n ; i++) { dest[i] = source[i]; }
 929       return dest;
 930     }
 931