PWGLF/QATasks/AliAnalysisTaskQAHighPtDeDx.cxx

   1  /*************************************************************************
   2 * Copyright(c) 1998-2008, 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 #include "AliAnalysisTaskQAHighPtDeDx.h"
  17
  18 // ROOT includes
  19 #include <TList.h>
  20 #include <TTree.h>
  21 #include <TMath.h>
  22 #include <TH1.h>
  23 #include <TF1.h>
  24 #include <TProfile.h>
  25 #include <TParticle.h>
  26 #include <TFile.h>
  27
  28 // AliRoot includes
  29 #include <AliAnalysisManager.h>
  30 #include <AliAnalysisFilter.h>
  31 #include <AliESDInputHandler.h>
  32 #include <AliESDEvent.h>
  33 #include <AliESDVertex.h>
  34 #include <AliLog.h>
  35 #include <AliExternalTrackParam.h>
  36 #include <AliESDtrackCuts.h>
  37 #include <AliESDVZERO.h>
  38 #include <AliAODVZERO.h>
  39
  40 #include <AliMCEventHandler.h>
  41 #include <AliMCEvent.h>
  42 #include <AliStack.h>
  43
  44 #include <TTreeStream.h>
  45
  46 #include <AliHeader.h>
  47 #include <AliGenPythiaEventHeader.h>
  48 #include <AliGenDPMjetEventHeader.h>
  49
  50 #include "AliCentrality.h"
  51 #include <AliESDv0.h>
  52 #include <AliKFVertex.h>
  53 #include <AliAODVertex.h>
  54
  55 #include <AliAODTrack.h>
  56 #include <AliAODPid.h>
  57 #include <AliAODMCHeader.h>
  58
  59
  60 // STL includes
  61 #include <iostream>
  62 using namespace std;
  63
  64
  65 //
  66 // Responsible:
  67 // Antonio Ortiz (Lund)
  68 // Peter Christiansen (Lund)
  69 //
  70
  71
  72
  73
  74
  75
  76
  77
  78 const Double_t AliAnalysisTaskQAHighPtDeDx::fgkClight = 2.99792458e-2;
  79 Float_t magf = -1;
  80 TF1* cutLow  = new TF1("StandardPhiCutLow",  "0.1/x/x+pi/18.0-0.025", 0, 50);
  81 TF1* cutHigh = new TF1("StandardPhiCutHigh", "0.12/x+pi/18.0+0.035", 0, 50);
  82 Double_t DeDxMIPMin  = 30;
  83 Double_t DeDxMIPMax  = 65;
  84 const Int_t nHists = 9;
  85 Float_t centralityGlobal = -10;
  86 Int_t etaLow[nHists]  = {-8, -8, -6, -4, -2, 0, 2, 4, 6};
  87 Int_t etaHigh[nHists] = { 8, -6, -4, -2,  0, 2, 4, 6, 8};
  88
  89 Int_t nDeltaPiBins   = 80;
  90 Double_t deltaPiLow  = 20;
  91 Double_t deltaPiHigh = 100;
  92 const Char_t *Pid[7]={"Ch","Pion","Kaon","Proton","Electron","Muon","Oher"};
  93 ClassImp(AliAnalysisTaskQAHighPtDeDx)
  94 //_____________________________________________________________________________
  95 //AliAnalysisTaskQAHighPtDeDx::AliAnalysisTaskQAHighPtDeDx(const char *name):
  96 AliAnalysisTaskQAHighPtDeDx::AliAnalysisTaskQAHighPtDeDx():
  97   AliAnalysisTaskSE(),
  98   fESD(0x0),
  99   fAOD(0x0),
 100   fMC(0x0),
 101   fMCStack(0x0),
 102   fMCArray(0x0),
 103   fTrackFilterGolden(0x0),
 104   fTrackFilterTPC(0x0),
 105   fCentEst("V0M"),
 106   fAnalysisType("ESD"),
 107   fAnalysisMC(kFALSE),
 108   fAnalysisPbPb(kFALSE),
 109   ftrigBit(0x0),
 110   fRandom(0x0),
 111   fPileUpRej(kFALSE),
 112   fVtxCut(10.0),
 113   fEtaCut(0.9),
 114   fMinCent(0.0),
 115   fMaxCent(100.0),
 116   fStoreMcIn(kFALSE),//
 117   fMcProcessType(-999),
 118   fTriggeredEventMB(-999),
 119   fVtxStatus(-999),
 120   fZvtx(-999),
 121   fZvtxMC(-999),
 122   fRun(-999),
 123   fEventId(-999),
 124   fListOfObjects(0),
 125   fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
 126   fn1(0x0),
 127   fcent(0x0),
 128   hMIPVsEta(0x0),
 129   pMIPVsEta(0x0),
 130   hMIPVsEtaV0s(0x0),
 131   pMIPVsEtaV0s(0x0),
 132   hPlateauVsEta(0x0),
 133   pPlateauVsEta(0x0),
 134   hPhi(0x0)
 135
 136
 137 {
 138   //default constructor
 139   for(Int_t i=0;i<9;++i){
 140
 141     hMIPVsNch[i]=0;//TH2D, MIP vs Nch for different eta intervals
 142     pMIPVsNch[i]=0;//TProfile, MIP vs Nch for different eta intervals
 143     hMIPVsPhi[i]=0;//TH2D, MIP vs phi for different eta intervals
 144     pMIPVsPhi[i]=0;//TProfile, MIP vs phi for different eta intervals
 145     hPlateauVsPhi[i]=0;//TH2D, dE/dx vs Phi, electrons 0.4<p<0.6 GeV/c
 146     pPlateauVsPhi[i]=0;//TProfile, dE/dx vs Phi, electrons 0.4<p<0.6 GeV/c
 147     histPiV0[i]=0;//TH2D, dE/dx vs p, pi id by V0s
 148     histpPiV0[i]=0;//TH1D, pi id by V0s
 149     histPV0[i]=0;// TH2D, dE/dx vs p, p id by V0s
 150     histpPV0[i]=0;// TH1D, p id by V0s
 151     histAllCh[i]=0;//TH2D, dE/dx vs p for all charged particles
 152     histPiTof[i]=0;//TH2D, dE/dx vs p for a "clean" sample of pions, beta>1
 153     histpPiTof[i]=0;//TH1D, for a "clean" sample of pions, beta>1
 154     histEV0[i]=0;
 155
 156     for(Int_t pid=0;pid<7;++pid){
 157       hMcIn[pid][i]=0;
 158       hMcOut[pid][i]=0;
 159     }
 160
 161   }
 162
 163
 164
 165 }
 166
 167
 168 AliAnalysisTaskQAHighPtDeDx::AliAnalysisTaskQAHighPtDeDx(const char *name):
 169   AliAnalysisTaskSE(name),
 170   fESD(0x0),
 171   fAOD(0x0),
 172   fMC(0x0),
 173   fMCStack(0x0),
 174   fMCArray(0x0),
 175   fTrackFilterGolden(0x0),
 176   fTrackFilterTPC(0x0),
 177   fCentEst("V0M"),
 178   fAnalysisType("ESD"),
 179   fAnalysisMC(kFALSE),
 180   fAnalysisPbPb(kFALSE),
 181   ftrigBit(0x0),
 182   fRandom(0x0),
 183   fPileUpRej(kFALSE),
 184   fVtxCut(10.0),
 185   fEtaCut(0.9),
 186   fMinCent(0.0),
 187   fMaxCent(100.0),
 188   fStoreMcIn(kFALSE),//
 189   fMcProcessType(-999),
 190   fTriggeredEventMB(-999),
 191   fVtxStatus(-999),
 192   fZvtx(-999),
 193   fZvtxMC(-999),
 194   fRun(-999),
 195   fEventId(-999),
 196   fListOfObjects(0),
 197   fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
 198   fn1(0x0),
 199   fcent(0x0),
 200   hMIPVsEta(0x0),
 201   pMIPVsEta(0x0),
 202   hMIPVsEtaV0s(0x0),
 203   pMIPVsEtaV0s(0x0),
 204   hPlateauVsEta(0x0),
 205   pPlateauVsEta(0x0),
 206   hPhi(0x0)
 207
 208
 209 {
 210   // Default constructor (should not be used)
 211   for(Int_t i=0;i<9;++i){
 212
 213     hMIPVsNch[i]=0;//TH2D, MIP vs Nch for different eta intervals
 214     pMIPVsNch[i]=0;//TProfile, MIP vs Nch for different eta intervals
 215     hMIPVsPhi[i]=0;//TH2D, MIP vs phi for different eta intervals
 216     pMIPVsPhi[i]=0;//TProfile, MIP vs phi for different eta intervals
 217     hPlateauVsPhi[i]=0;//TH2D, dE/dx vs Phi, electrons 0.4<p<0.6 GeV/c
 218     pPlateauVsPhi[i]=0;//TProfile, dE/dx vs Phi, electrons 0.4<p<0.6 GeV/c
 219     histPiV0[i]=0;//TH2D, dE/dx vs p, pi id by V0s
 220     histpPiV0[i]=0;//TH1D, pi id by V0s
 221     histPV0[i]=0;// TH2D, dE/dx vs p, p id by V0s
 222     histpPV0[i]=0;// TH1D, p id by V0s
 223     histAllCh[i]=0;//TH2D, dE/dx vs p for all charged particles
 224     histPiTof[i]=0;//TH2D, dE/dx vs p for a "clean" sample of pions, beta>1
 225     histpPiTof[i]=0;//TH1D, for a "clean" sample of pions, beta>1
 226     histEV0[i]=0;
 227
 228     for(Int_t pid=0;pid<7;++pid){
 229       hMcIn[pid][i]=0;
 230       hMcOut[pid][i]=0;
 231     }
 232
 233   }
 234   DefineOutput(1, TList::Class());//esto es nuevo
 235 }
 236
 237
 238
 239
 240 AliAnalysisTaskQAHighPtDeDx::~AliAnalysisTaskQAHighPtDeDx() {
 241   //
 242   // Destructor
 243   //
 244
 245 }
 246
 247
 248
 249
 250
 251
 252
 253
 254 //______________________________________________________________________________
 255 void AliAnalysisTaskQAHighPtDeDx::UserCreateOutputObjects()
 256 {
 257   // This method is called once per worker node
 258   // Here we define the output: histograms and debug tree if requested
 259   // We also create the random generator here so it might get different seeds...
 260   fRandom = new TRandom(0); // 0 means random seed
 261
 262
 263
 264   //OpenFile(1);
 265   fListOfObjects = new TList();
 266   fListOfObjects->SetOwner();
 267
 268
 269
 270
 271   //
 272   // Histograms
 273   //
 274   fEvents = new TH1I("fEvents","Number of analyzed events; Events; Counts", 3, 0, 3);
 275   fListOfObjects->Add(fEvents);
 276
 277   fn1=new TH1F("fn1","fn1",11,-1,10);
 278   fListOfObjects->Add(fn1);
 279
 280   fcent=new TH1F("fcent","fcent",104,-2,102);
 281   fListOfObjects->Add(fcent);
 282
 283   fVtx = new TH1I("fVtx","Vtx info (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
 284   fListOfObjects->Add(fVtx);
 285
 286   fVtxBeforeCuts = new TH1F("fVtxBeforeCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
 287   fListOfObjects->Add(fVtxBeforeCuts);
 288
 289   fVtxAfterCuts = new TH1F("fVtxAfterCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
 290   fListOfObjects->Add(fVtxAfterCuts);
 291
 292
 293   const Int_t nPtBinsV0s = 25;
 294   Double_t ptBinsV0s[nPtBinsV0s+1] = { 0.0 , 0.1 , 0.2 , 0.3 , 0.4 , 0.5 , 0.6 , 0.7 , 0.8 , 0.9 , 1.0 ,
 295                                        1.2 , 1.4 , 1.6 , 1.8 , 2.0 , 2.5 , 3.0 , 3.5 , 4.0 , 5.0 , 7.0 ,
 296                                        9.0 , 12.0, 15.0, 20.0 };
 297
 298
 299
 300
 301   const Char_t* ending[nHists] = {"", "86", "64", "42", "20", "02", "24", "46", "68"};
 302
 303   const Char_t* LatexEta[nHists] = {
 304     "|#eta|<0.8", "-0.8<#eta<-0.6", "-0.6<#eta<-0.4", "-0.4<#eta<-0.2", "-0.2<#eta<0",
 305     "0<#eta<0.2", "0.2<#eta<0.4", "0.4<#eta<0.6", "0.6<#eta<0.8"
 306   };
 307   hPhi = new TH2D("histPhi", "pt; #phi'", nPtBinsV0s, ptBinsV0s, 90, -0.05, 0.4);
 308   //dE/dx vs phi, pions at the MIP
 309   fListOfObjects->Add(hPhi);
 310
 311
 312
 313
 314   Int_t nPhiBins = 36;
 315
 316   for(Int_t i = 0; i < nHists; i++) {
 317
 318
 319     hMIPVsPhi[i]      = new TH2D(Form("hMIPVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx MIP",LatexEta[i]), nPhiBins, 0, 2*TMath::Pi(),
 320                                  DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
 321     hMIPVsPhi[i]->Sumw2();
 322
 323     pMIPVsPhi[i]     = new TProfile(Form("pMIPVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx MIP",LatexEta[i]),  nPhiBins, 0, 2*TMath::Pi(),
 324                                     DeDxMIPMin, DeDxMIPMax);
 325     pMIPVsPhi[i]->SetMarkerStyle(20);
 326     pMIPVsPhi[i]->SetMarkerColor(1);
 327     pMIPVsPhi[i]->SetLineColor(1);
 328     pMIPVsPhi[i]->Sumw2();
 329
 330     fListOfObjects->Add(hMIPVsPhi[i]);
 331     fListOfObjects->Add(pMIPVsPhi[i]);
 332
 333     hPlateauVsPhi[i]  = new TH2D(Form("hPlateauVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx Plateau",LatexEta[i]),  nPhiBins, 0, 2*TMath::Pi(),
 334                                   95-DeDxMIPMax, DeDxMIPMax, 95);
 335     hPlateauVsPhi[i]->Sumw2();
 336
 337     pPlateauVsPhi[i] = new TProfile(Form("pPlateauVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx Plateau",LatexEta[i]), nPhiBins, 0, 2*TMath::Pi(),
 338                                      DeDxMIPMax, 95);
 339     pPlateauVsPhi[i]->SetMarkerStyle(20);
 340     pPlateauVsPhi[i]->SetMarkerColor(1);
 341     pPlateauVsPhi[i]->SetLineColor(1);
 342     pPlateauVsPhi[i]->Sumw2();
 343
 344     fListOfObjects->Add(hPlateauVsPhi[i]);
 345     fListOfObjects->Add(pPlateauVsPhi[i]);
 346
 347
 348     hMIPVsNch[i]      = new TH2D(Form("hMIPVsNch%s", ending[i]), Form("%s; TPC track mult. |#eta|<0.8; dE/dx MIP",LatexEta[i]), 400, 1, 2001,
 349                                   DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
 350     hMIPVsNch[i]->Sumw2();
 351
 352     pMIPVsNch[i] = new TProfile(Form("pMIPVsNch%s", ending[i]), Form("%s; TPC track mult. |#eta|<0.8; dE/dx MIP",LatexEta[i]), 400, 1, 2001, DeDxMIPMin, DeDxMIPMax);
 353     pMIPVsNch[i]->SetMarkerStyle(20);
 354     pMIPVsNch[i]->SetMarkerColor(1);
 355     pMIPVsNch[i]->SetLineColor(1);
 356     pMIPVsNch[i]->Sumw2();
 357
 358     fListOfObjects->Add(hMIPVsNch[i]);
 359     fListOfObjects->Add(pMIPVsNch[i]);
 360
 361     //two dimmesional distributions dE/dx vs p for secondary pions
 362     histPiV0[i]  = new TH2D(Form("histPiV0%s", ending[i]), "Pions id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
 363     histPiV0[i]->Sumw2();
 364     fListOfObjects->Add(histPiV0[i]);
 365
 366     histpPiV0[i]  = new TH1D(Form("histPiV01D%s", ending[i]), "Pions id by V0; #it{p} (GeV/#it{c}); counts", 200, 0, 20);
 367     histpPiV0[i]->Sumw2();
 368     fListOfObjects->Add(histpPiV0[i]);
 369
 370     //two dimmesional distributions dE/dx vs p for secondary protons
 371     histPV0[i]   = new TH2D(Form("histPV0%s", ending[i]), "Protons id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
 372     histPV0[i]->Sumw2();
 373     fListOfObjects->Add(histPV0[i]);
 374
 375     histpPV0[i]  = new TH1D(Form("histPV01D%s", ending[i]), "Protons id by V0; #it{p} (GeV/#it{c}); counts", 200, 0, 20);
 376     histpPV0[i]->Sumw2();
 377     fListOfObjects->Add(histpPV0[i]);
 378
 379     //two dimmesional distributions dE/dx vs p for primary pions
 380     histPiTof[i] = new TH2D(Form("histPiTof%s", ending[i]), "all charged", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
 381     histPiTof[i]->Sumw2();
 382
 383     histpPiTof[i]  = new TH1D(Form("histPiTof1D%s", ending[i]), "Protons id by V0; #it{p} (GeV/#it{c}); counts", 200, 0, 20);
 384     histpPiTof[i]->Sumw2();
 385     fListOfObjects->Add(histpPiTof[i]);
 386
 387
 388     histAllCh[i] = new TH2D(Form("histAllCh%s", ending[i]), "Pions id by TOF", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
 389     histAllCh[i]->Sumw2();
 390
 391     fListOfObjects->Add(histPiTof[i]);
 392     fListOfObjects->Add(histAllCh[i]);
 393
 394     histEV0[i]   = new TH2D(Form("histEV0%s", ending[i]), "Electrons id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
 395     histEV0[i]->Sumw2();
 396     fListOfObjects->Add(histEV0[i]);
 397
 398
 399
 400   }
 401
 402
 403   hMIPVsEta = new TH2D("hMIPVsEta","; #eta; dE/dx_{MIP, primary tracks}",16,-0.8,0.8,DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
 404   pMIPVsEta = new TProfile("pMIPVsEta","; #eta; #LT dE/dx #GT_{MIP, primary tracks}",16,-0.8,0.8, DeDxMIPMin, DeDxMIPMax);
 405   hMIPVsEtaV0s = new TH2D("hMIPVsEtaV0s","; #eta; dE/dx_{MIP, secondary tracks}",16,-0.8,0.8,DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
 406   pMIPVsEtaV0s = new TProfile("pMIPVsEtaV0s","; #eta; #LT dE/dx #GT_{MIP, secondary tracks}",16,-0.8,0.8,DeDxMIPMin, DeDxMIPMax);
 407
 408   hPlateauVsEta = new TH2D("hPlateauVsEta","; #eta; dE/dx_{Plateau, primary tracks}",16,-0.8,0.8,95-DeDxMIPMax, DeDxMIPMax, 95);
 409   pPlateauVsEta = new TProfile("pPlateauVsEta","; #eta; #LT dE/dx #GT_{Plateau, primary tracks}",16,-0.8,0.8, DeDxMIPMax, 95);
 410
 411   fListOfObjects->Add(hMIPVsEta);
 412   fListOfObjects->Add(pMIPVsEta);
 413   fListOfObjects->Add(hMIPVsEtaV0s);
 414   fListOfObjects->Add(pMIPVsEtaV0s);
 415   fListOfObjects->Add(hPlateauVsEta);
 416   fListOfObjects->Add(pPlateauVsEta);
 417
 418
 419
 420
 421
 422
 423   if (fAnalysisMC) {
 424     for(Int_t i = 0; i < nHists; i++) {
 425       for(Int_t pid = 0; pid < 7; pid++) {
 426
 427         hMcIn[pid][i] = new TH1D(Form("hIn%s%s", Pid[pid],ending[i]), Form("MC in (pid %s)", Pid[pid]),
 428                                  nPtBinsV0s, ptBinsV0s);
 429         fListOfObjects->Add(hMcIn[pid][i]);
 430
 431         hMcOut[pid][i] = new TH1D(Form("hMcOut%s%s", Pid[pid],ending[i]), Form("MC out (pid %s)", Pid[pid]),
 432                                   nPtBinsV0s, ptBinsV0s);
 433         fListOfObjects->Add(hMcOut[pid][i]);
 434
 435
 436       }
 437     }
 438
 439     fVtxMC = new TH1F("fVtxMC","mc vtx", 120, -30, 30);
 440     fListOfObjects->Add(fVtxMC);
 441
 442   }
 443
 444   // Post output data.
 445   PostData(1, fListOfObjects);
 446 }
 447
 448 //______________________________________________________________________________
 449 void AliAnalysisTaskQAHighPtDeDx::UserExec(Option_t *)
 450 {
 451   // Main loop
 452
 453   //
 454   // First we make sure that we have valid input(s)!
 455   //
 456
 457
 458
 459   AliVEvent *event = InputEvent();
 460   if (!event) {
 461      Error("UserExec", "Could not retrieve event");
 462      return;
 463   }
 464
 465
 466
 467   if (fAnalysisType == "ESD"){
 468     fESD = dynamic_cast<AliESDEvent*>(event);
 469     if(!fESD){
 470       Printf("%s:%d ESDEvent not found in Input Manager",(char*)__FILE__,__LINE__);
 471       this->Dump();
 472       return;
 473     }
 474   } else {
 475     fAOD = dynamic_cast<AliAODEvent*>(event);
 476     if(!fAOD){
 477       Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
 478       this->Dump();
 479       return;
 480     }
 481   }
 482
 483
 484
 485   if (fAnalysisMC) {
 486
 487     if (fAnalysisType == "ESD"){
 488       fMC = dynamic_cast<AliMCEvent*>(MCEvent());
 489       if(!fMC){
 490         Printf("%s:%d MCEvent not found in Input Manager",(char*)__FILE__,__LINE__);
 491         this->Dump();
 492         return;
 493       }
 494
 495       fMCStack = fMC->Stack();
 496
 497       if(!fMCStack){
 498         Printf("%s:%d MCStack not found in Input Manager",(char*)__FILE__,__LINE__);
 499         this->Dump();
 500         return;
 501       }
 502     } else { // AOD
 503
 504       fMC = dynamic_cast<AliMCEvent*>(MCEvent());
 505       if(fMC)
 506         fMC->Dump();
 507
 508       fMCArray = (TClonesArray*)fAOD->FindListObject("mcparticles");
 509       if(!fMCArray){
 510         Printf("%s:%d AOD MC array not found in Input Manager",(char*)__FILE__,__LINE__);
 511         this->Dump();
 512         return;
 513       }
 514     }
 515   }
 516
 517
 518   // Get trigger decision
 519   fTriggeredEventMB = 0; //init
 520
 521   fn1->Fill(0);
 522
 523   if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 524      ->IsEventSelected() & ftrigBit ){
 525     fTriggeredEventMB = 1;  //event triggered as minimum bias
 526   }
 527
 528   // Get process type for MC
 529   fMcProcessType = 0; // -1=invalid, 0=data, 1=ND, 2=SD, 3=DD
 530
 531   // real data that are not triggered we skip
 532   if(!fAnalysisMC && !fTriggeredEventMB)
 533     return;
 534
 535   fn1->Fill(1);
 536
 537
 538   if (fAnalysisMC) {
 539
 540
 541
 542     if (fAnalysisType == "ESD"){
 543
 544
 545
 546       AliHeader* headerMC = fMC->Header();
 547       if (headerMC) {
 548
 549         AliGenEventHeader* genHeader = headerMC->GenEventHeader();
 550         TArrayF vtxMC(3); // primary vertex  MC
 551         vtxMC[0]=9999; vtxMC[1]=9999;  vtxMC[2]=9999; //initialize with dummy
 552         if (genHeader) {
 553           genHeader->PrimaryVertex(vtxMC);
 554         }
 555         fZvtxMC = vtxMC[2];
 556
 557         // PYTHIA:
 558         AliGenPythiaEventHeader* pythiaGenHeader =
 559           dynamic_cast<AliGenPythiaEventHeader*>(headerMC->GenEventHeader());
 560         if (pythiaGenHeader) {  //works only for pythia
 561           fMcProcessType =  GetPythiaEventProcessType(pythiaGenHeader->ProcessType());
 562         }
 563         // PHOJET:
 564         AliGenDPMjetEventHeader* dpmJetGenHeader =
 565           dynamic_cast<AliGenDPMjetEventHeader*>(headerMC->GenEventHeader());
 566         if (dpmJetGenHeader) {
 567           fMcProcessType = GetDPMjetEventProcessType(dpmJetGenHeader->ProcessType());
 568         }
 569       }
 570     } else { // AOD
 571
 572
 573
 574       AliAODMCHeader* mcHeader = dynamic_cast<AliAODMCHeader*>(fAOD->FindListObject("mcHeader"));
 575
 576
 577       if(mcHeader) {
 578         fZvtxMC = mcHeader->GetVtxZ();
 579
 580
 581
 582         if(strstr(mcHeader->GetGeneratorName(), "Pythia")) {
 583           fMcProcessType =  GetPythiaEventProcessType(mcHeader->GetEventType());
 584         } else {
 585           fMcProcessType =  GetDPMjetEventProcessType(mcHeader->GetEventType());
 586         }
 587       }
 588     }
 589
 590
 591   }
 592
 593
 594
 595   if (fAnalysisType == "ESD"){
 596
 597     const AliESDVertex *vtxESD = fESD->GetPrimaryVertexTracks();
 598     if(vtxESD->GetNContributors()<1) {
 599       // SPD vertex
 600       vtxESD = fESD->GetPrimaryVertexSPD();
 601       /* quality checks on SPD-vertex */
 602       TString vertexType = vtxESD->GetTitle();
 603       if (vertexType.Contains("vertexer: Z") && (vtxESD->GetDispersion() > 0.04 || vtxESD->GetZRes() > 0.25))
 604         fZvtx  = -1599; //vertex = 0x0; //
 605       else if (vtxESD->GetNContributors()<1)
 606         fZvtx  = -999; //vertex = 0x0; //
 607       else
 608         fZvtx = vtxESD->GetZ();
 609     }
 610     else
 611       fZvtx = vtxESD->GetZ();
 612
 613   }
 614   else // AOD
 615     fZvtx = GetVertex(fAOD);
 616
 617   fVtxBeforeCuts->Fill(fZvtx);
 618
 619   //cut on the z position of vertex
 620   if (TMath::Abs(fZvtx) > fVtxCut) {
 621     return;
 622   }
 623   fn1->Fill(2);
 624
 625
 626
 627
 628
 629
 630   Float_t centrality = -10;
 631
 632   // only analyze triggered events
 633   if(fTriggeredEventMB) {
 634
 635     if (fAnalysisType == "ESD"){
 636       if(fAnalysisPbPb){
 637         AliCentrality *centObject = fESD->GetCentrality();
 638         centrality = centObject->GetCentralityPercentile(fCentEst);
 639         centralityGlobal = centrality;
 640         if((centrality>fMaxCent)||(centrality<fMinCent))return;
 641       }
 642       fcent->Fill(centrality);
 643       fn1->Fill(3);
 644       if(fAnalysisMC){
 645         if(TMath::Abs(fZvtxMC)<fVtxCut){
 646           ProcessMCTruthESD();
 647           fVtxMC->Fill(fZvtxMC);
 648         }
 649       }
 650       AnalyzeESD(fESD);
 651     } else { // AOD
 652       if(fAnalysisPbPb){
 653         AliCentrality *centObject = fAOD->GetCentrality();
 654         if(centObject){
 655           centrality = centObject->GetCentralityPercentile(fCentEst);
 656         }
 657         //cout<<"centrality="<<centrality<<endl;
 658         if((centrality>fMaxCent)||(centrality<fMinCent))return;
 659       }
 660       fcent->Fill(centrality);
 661       fn1->Fill(3);
 662       if(fAnalysisMC){
 663         if(TMath::Abs(fZvtxMC)<fVtxCut){
 664
 665           ProcessMCTruthAOD();
 666           fVtxMC->Fill(fZvtxMC);
 667         }
 668       }
 669       AnalyzeAOD(fAOD);
 670     }
 671   }
 672
 673   fVtxAfterCuts->Fill(fZvtx);
 674
 675
 676
 677
 678   // Post output data.
 679   PostData(1, fListOfObjects);
 680 }
 681
 682 //________________________________________________________________________
 683 void AliAnalysisTaskQAHighPtDeDx::AnalyzeESD(AliESDEvent* esdEvent)
 684 {
 685   fRun  = esdEvent->GetRunNumber();
 686   fEventId = 0;
 687   if(esdEvent->GetHeader())
 688     fEventId = GetEventIdAsLong(esdEvent->GetHeader());
 689
 690   cout << "centrality=" << centralityGlobal << endl;
 691
 692
 693   Bool_t isPileup = esdEvent->IsPileupFromSPD();
 694   if(fPileUpRej)
 695     if(isPileup)
 696       return;
 697   fn1->Fill(4);
 698
 699
 700   //  Int_t     event     = esdEvent->GetEventNumberInFile();
 701   //UInt_t    time      = esdEvent->GetTimeStamp();
 702   //  ULong64_t trigger   = esdEvent->GetTriggerMask();
 703   magf      = esdEvent->GetMagneticField();
 704
 705
 706
 707
 708
 709   if(fTriggeredEventMB) {// Only MC case can we have not triggered events
 710
 711     // accepted event
 712     fEvents->Fill(0);
 713
 714
 715     //Change, 10/04/13. Now accept all events to do a correct normalization
 716     //if(fVtxStatus!=1) return; // accepted vertex
 717     //Int_t nESDTracks = esdEvent->GetNumberOfTracks();
 718
 719     ProduceArrayTrksESD( esdEvent );//produce array with global track parameters
 720     ProduceArrayV0ESD( esdEvent );//v0's
 721
 722
 723     fEvents->Fill(1);
 724
 725
 726
 727
 728   } // end if triggered
 729
 730
 731 }
 732
 733 //________________________________________________________________________
 734 void AliAnalysisTaskQAHighPtDeDx::AnalyzeAOD(AliAODEvent* aodEvent)
 735 {
 736   fRun  = aodEvent->GetRunNumber();
 737   fEventId = 0;
 738   if(aodEvent->GetHeader())
 739     fEventId = GetEventIdAsLong(aodEvent->GetHeader());
 740
 741   //UInt_t    time      = 0; // Missing AOD info? aodEvent->GetTimeStamp();
 742   magf      = aodEvent->GetMagneticField();
 743
 744   //Int_t     trackmult = 0; // no pt cuts
 745   //Int_t     nadded    = 0;
 746
 747   Bool_t isPileup = aodEvent->IsPileupFromSPD();
 748   if(fPileUpRej)
 749     if(isPileup)
 750       return;
 751   fn1->Fill(4);
 752
 753
 754
 755   if(fTriggeredEventMB) {// Only MC case can we have not triggered events
 756
 757     // accepted event
 758     fEvents->Fill(0);
 759
 760     //if(fVtxStatus!=1) return; // accepted vertex
 761     //Int_t nAODTracks = aodEvent->GetNumberOfTracks();
 762
 763     ProduceArrayTrksAOD( aodEvent );
 764     ProduceArrayV0AOD( aodEvent );
 765
 766     fEvents->Fill(1);
 767
 768
 769
 770
 771   } // end if triggered
 772
 773 }
 774
 775 //_____________________________________________________________________________
 776 Float_t AliAnalysisTaskQAHighPtDeDx::GetVertex(const AliVEvent* event) const
 777 {
 778   Float_t zvtx = -999;
 779
 780   const AliVVertex* primaryVertex = event->GetPrimaryVertex();
 781
 782   if(primaryVertex->GetNContributors()>0)
 783     zvtx = primaryVertex->GetZ();
 784
 785   return zvtx;
 786 }
 787
 788 //_____________________________________________________________________________
 789 Short_t AliAnalysisTaskQAHighPtDeDx::GetPidCode(Int_t pdgCode) const
 790 {
 791   // return our internal code for pions, kaons, and protons
 792
 793   Short_t pidCode = 6;
 794
 795   switch (TMath::Abs(pdgCode)) {
 796   case 211:
 797     pidCode = 1; // pion
 798     break;
 799   case 321:
 800     pidCode = 2; // kaon
 801     break;
 802   case 2212:
 803     pidCode = 3; // proton
 804     break;
 805   case 11:
 806     pidCode = 4; // electron
 807     break;
 808   case 13:
 809     pidCode = 5; // muon
 810     break;
 811   default:
 812     pidCode = 6;  // something else?
 813   };
 814
 815   return pidCode;
 816 }
 817
 818 //_____________________________________________________________________________
 819 void AliAnalysisTaskQAHighPtDeDx::ProcessMCTruthESD()
 820 {
 821   // Fill the special MC histogram with the MC truth info
 822
 823   const Int_t nTracksMC = fMCStack->GetNtrack();
 824
 825   for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
 826
 827     //Cuts
 828     if(!(fMCStack->IsPhysicalPrimary(iTracks)))
 829       continue;
 830
 831     TParticle* trackMC = fMCStack->Particle(iTracks);
 832
 833     TParticlePDG* pdgPart = trackMC ->GetPDG();
 834     Double_t chargeMC = pdgPart->Charge();
 835
 836     if(chargeMC==0)
 837       continue;
 838
 839     if (TMath::Abs(trackMC->Eta()) > fEtaCut )
 840       continue;
 841
 842     Double_t etaMC = trackMC->Eta();
 843     Int_t pdgCode = trackMC->GetPdgCode();
 844     Short_t pidCodeMC = 0;
 845     pidCodeMC = GetPidCode(pdgCode);
 846
 847
 848     for(Int_t nh = 0; nh < 9; nh++) {
 849
 850       if( etaMC > etaHigh[nh]/10.0 || etaMC < etaLow[nh]/10.0 )
 851         continue;
 852
 853       hMcIn[0][nh]->Fill(trackMC->Pt());
 854       hMcIn[pidCodeMC][nh]->Fill(trackMC->Pt());
 855
 856
 857     }
 858
 859   }//MC track loop
 860
 861
 862
 863 }
 864
 865 //_____________________________________________________________________________
 866 void AliAnalysisTaskQAHighPtDeDx::ProcessMCTruthAOD()
 867 {
 868   // Fill the special MC histogram with the MC truth info
 869
 870
 871   const Int_t nTracksMC = fMCArray->GetEntriesFast();
 872
 873   for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
 874
 875     AliAODMCParticle* trackMC = dynamic_cast<AliAODMCParticle*>(fMCArray->At(iTracks));
 876
 877     if(!trackMC){
 878       AliError("Cannot get MC particle");
 879       continue;
 880     }
 881
 882     //Cuts
 883     if(!(trackMC->IsPhysicalPrimary()))
 884       continue;
 885
 886
 887     Double_t chargeMC = trackMC->Charge();
 888     if(chargeMC==0)
 889       continue;
 890
 891
 892     if (TMath::Abs(trackMC->Eta()) > fEtaCut )
 893       continue;
 894
 895     Double_t etaMC = trackMC->Eta();
 896     Int_t pdgCode = trackMC->GetPdgCode();
 897     Short_t pidCodeMC = 0;
 898     pidCodeMC = GetPidCode(pdgCode);
 899
 900     //cout<<"pidcode="<<pidCodeMC<<endl;
 901     for(Int_t nh = 0; nh < 9; nh++) {
 902
 903       if( etaMC > etaHigh[nh]/10.0 || etaMC < etaLow[nh]/10.0 )
 904         continue;
 905
 906       hMcIn[0][nh]->Fill(trackMC->Pt());
 907       hMcIn[pidCodeMC][nh]->Fill(trackMC->Pt());
 908
 909
 910     }
 911
 912   }//MC track loop
 913
 914
 915 }
 916
 917
 918 //_____________________________________________________________________________
 919 Short_t AliAnalysisTaskQAHighPtDeDx::GetPythiaEventProcessType(Int_t pythiaType) {
 920   //
 921   // Get the process type of the event.  PYTHIA
 922   //
 923   // source PWG0   dNdpt
 924
 925   Short_t globalType = -1; //init
 926
 927   if(pythiaType==92||pythiaType==93){
 928     globalType = 2; //single diffractive
 929   }
 930   else if(pythiaType==94){
 931     globalType = 3; //double diffractive
 932   }
 933   //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??
 934   else {
 935     globalType = 1;  //non diffractive
 936   }
 937   return globalType;
 938 }
 939
 940 //_____________________________________________________________________________
 941 Short_t AliAnalysisTaskQAHighPtDeDx::GetDPMjetEventProcessType(Int_t dpmJetType) {
 942   //
 943   // get the process type of the event.  PHOJET
 944   //
 945   //source PWG0   dNdpt
 946   // can only read pythia headers, either directly or from cocktalil header
 947   Short_t globalType = -1;
 948
 949   if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction
 950     globalType = 1;
 951   }
 952   else if (dpmJetType==5 || dpmJetType==6) {
 953     globalType = 2;
 954   }
 955   else if (dpmJetType==7) {
 956     globalType = 3;
 957   }
 958   return globalType;
 959 }
 960
 961 //_____________________________________________________________________________
 962 ULong64_t AliAnalysisTaskQAHighPtDeDx::GetEventIdAsLong(AliVHeader* header) const
 963 {
 964   // To have a unique id for each event in a run!
 965   // Modified from AliRawReader.h
 966   return ((ULong64_t)header->GetBunchCrossNumber()+
 967           (ULong64_t)header->GetOrbitNumber()*3564+
 968           (ULong64_t)header->GetPeriodNumber()*16777215*3564);
 969 }
 970
 971
 972 //____________________________________________________________________
 973 TParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMother(AliStack* stack, Int_t label)
 974 {
 975   //
 976   // Finds the first mother among the primary particles of the particle identified by <label>,
 977   // i.e. the primary that "caused" this particle
 978   //
 979   // Taken from AliPWG0Helper class
 980   //
 981
 982   Int_t motherLabel = FindPrimaryMotherLabel(stack, label);
 983   if (motherLabel < 0)
 984     return 0;
 985
 986   return stack->Particle(motherLabel);
 987 }
 988
 989 //____________________________________________________________________
 990 Int_t AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherLabel(AliStack* stack, Int_t label)
 991 {
 992   //
 993   // Finds the first mother among the primary particles of the particle identified by <label>,
 994   // i.e. the primary that "caused" this particle
 995   //
 996   // returns its label
 997   //
 998   // Taken from AliPWG0Helper class
 999   //
1000   const Int_t nPrim  = stack->GetNprimary();
1001
1002   while (label >= nPrim) {
1003
1004     //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
1005
1006     TParticle* particle = stack->Particle(label);
1007     if (!particle) {
1008
1009       AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
1010       return -1;
1011     }
1012
1013     // find mother
1014     if (particle->GetMother(0) < 0) {
1015
1016       AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
1017       return -1;
1018     }
1019
1020     label = particle->GetMother(0);
1021   }
1022
1023   return label;
1024 }
1025
1026 //____________________________________________________________________
1027 AliAODMCParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherAOD(AliAODMCParticle* startParticle)
1028 {
1029   //
1030   // Finds the first mother among the primary particles of the particle identified by <label>,
1031   // i.e. the primary that "caused" this particle
1032   //
1033   // Taken from AliPWG0Helper class
1034   //
1035
1036   AliAODMCParticle* mcPart = startParticle;
1037
1038   while (mcPart)
1039     {
1040
1041       if(mcPart->IsPrimary())
1042         return mcPart;
1043
1044       Int_t mother = mcPart->GetMother();
1045
1046       mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));
1047     }
1048
1049   return 0;
1050 }
1051
1052
1053 //V0______________________________________
1054 //____________________________________________________________________
1055 TParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherV0(AliStack* stack, Int_t label)
1056 {
1057   //
1058   // Finds the first mother among the primary particles of the particle identified by <label>,
1059   // i.e. the primary that "caused" this particle
1060   //
1061   // Taken from AliPWG0Helper class
1062   //
1063
1064   Int_t nSteps = 0;
1065
1066   Int_t motherLabel = FindPrimaryMotherLabelV0(stack, label, nSteps);
1067   if (motherLabel < 0)
1068     return 0;
1069
1070   return stack->Particle(motherLabel);
1071 }
1072
1073 //____________________________________________________________________
1074 Int_t AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherLabelV0(AliStack* stack, Int_t label, Int_t& nSteps)
1075 {
1076   //
1077   // Finds the first mother among the primary particles of the particle identified by <label>,
1078   // i.e. the primary that "caused" this particle
1079   //
1080   // returns its label
1081   //
1082   // Taken from AliPWG0Helper class
1083   //
1084   nSteps = 0;
1085   const Int_t nPrim  = stack->GetNprimary();
1086
1087   while (label >= nPrim) {
1088
1089     //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
1090
1091     nSteps++; // 1 level down
1092
1093     TParticle* particle = stack->Particle(label);
1094     if (!particle) {
1095
1096       AliDebugGeneral("FindPrimaryMotherLabelV0", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
1097       return -1;
1098     }
1099
1100     // find mother
1101     if (particle->GetMother(0) < 0) {
1102
1103       AliDebugGeneral("FindPrimaryMotherLabelV0", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
1104       return -1;
1105     }
1106
1107     label = particle->GetMother(0);
1108   }
1109
1110   return label;
1111 }
1112
1113 //____________________________________________________________________
1114 AliAODMCParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherAODV0(AliAODMCParticle* startParticle, Int_t& nSteps)
1115 {
1116   //
1117   // Finds the first mother among the primary particles of the particle identified by <label>,
1118   // i.e. the primary that "caused" this particle
1119   //
1120   // Taken from AliPWG0Helper class
1121   //
1122
1123   nSteps = 0;
1124
1125   AliAODMCParticle* mcPart = startParticle;
1126
1127   while (mcPart)
1128     {
1129
1130       if(mcPart->IsPrimary())
1131         return mcPart;
1132
1133       Int_t mother = mcPart->GetMother();
1134
1135       mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));
1136       nSteps++; // 1 level down
1137     }
1138
1139   return 0;
1140 }
1141
1142
1143
1144 //__________________________________________________________________
1145 void AliAnalysisTaskQAHighPtDeDx::ProduceArrayTrksESD( AliESDEvent *ESDevent ){
1146
1147   const Int_t nESDTracks = ESDevent->GetNumberOfTracks();
1148   //Int_t trackmult=0;
1149
1150
1151   Int_t multTPC = 0;
1152
1153   //get multiplicity tpc only track cuts
1154   for(Int_t iT = 0; iT < nESDTracks; iT++) {
1155
1156     AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
1157
1158
1159     if(TMath::Abs(esdTrack->Eta()) > fEtaCut)
1160       continue;
1161
1162     //only golden track cuts
1163     UInt_t selectDebug = 0;
1164     if (fTrackFilterTPC) {
1165       selectDebug = fTrackFilterTPC->IsSelected(esdTrack);
1166       if (!selectDebug) {
1167         //cout<<"this is not a golden track"<<endl;
1168         continue;
1169       }
1170     }
1171
1172     multTPC++;
1173
1174   }
1175
1176
1177
1178   for(Int_t iT = 0; iT < nESDTracks; iT++) {
1179
1180     AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
1181
1182
1183     if(TMath::Abs(esdTrack->Eta()) > fEtaCut)
1184       continue;
1185
1186     //only golden track cuts
1187     UInt_t selectDebug = 0;
1188     if (fTrackFilterGolden) {
1189       selectDebug = fTrackFilterGolden->IsSelected(esdTrack);
1190       if (!selectDebug) {
1191         //cout<<"this is not a golden track"<<endl;
1192         continue;
1193       }
1194     }
1195
1196
1197     //Int_t sharedtpcclusters = esdTrack->GetTPCnclsS();
1198     Short_t ncl     = esdTrack->GetTPCsignalN();
1199
1200
1201     if(ncl<70)
1202       continue;
1203     Double_t eta  = esdTrack->Eta();
1204     Double_t phi  = esdTrack->Phi();
1205     Double_t momentum = esdTrack->P();
1206     Float_t  dedx    = esdTrack->GetTPCsignal();
1207
1208     //cout<<"magf="<<magf<<endl;
1209
1210     if(!PhiCut(esdTrack->Pt(), phi, esdTrack->Charge(), magf, cutLow, cutHigh))
1211       continue;
1212
1213
1214
1215
1216     //TOF
1217
1218     Bool_t IsTOFout=kFALSE;
1219     if ((esdTrack->GetStatus()&AliESDtrack::kTOFout)==0)
1220       IsTOFout=kTRUE;
1221     Float_t lengthtrack=esdTrack->GetIntegratedLength();
1222     Float_t timeTOF=esdTrack->GetTOFsignal();
1223     Double_t inttime[5]={0,0,0,0,0};
1224     esdTrack->GetIntegratedTimes(inttime);// Returns the array with integrated times for each particle hypothesis
1225     Float_t beta = -99;
1226     if ( !IsTOFout ){
1227       if ( ( lengthtrack != 0 ) && ( timeTOF != 0) )
1228         beta = inttime[0] / timeTOF;
1229     }
1230
1231     Short_t pidCode     = 0;
1232
1233     if(fAnalysisMC) {
1234
1235       const Int_t label = TMath::Abs(esdTrack->GetLabel());
1236       TParticle* mcTrack = fMCStack->Particle(label);
1237       if (mcTrack){
1238
1239         Int_t pdgCode = mcTrack->GetPdgCode();
1240         pidCode = GetPidCode(pdgCode);
1241
1242       }
1243
1244     }
1245
1246
1247     if( momentum <= 0.6 && momentum >= 0.4  ){//only p:0.4-0.6 GeV, pion MIP
1248
1249       if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
1250         if(momentum<0.6&&momentum>0.4){
1251           hMIPVsEta->Fill(eta,dedx);
1252           pMIPVsEta->Fill(eta,dedx);
1253         }
1254       }
1255       if( dedx > DeDxMIPMax+1 && dedx < 95 ){
1256         if(TMath::Abs(beta-1)<0.1){
1257           hPlateauVsEta->Fill(eta,dedx);
1258           pPlateauVsEta->Fill(eta,dedx);
1259         }
1260       }
1261     }
1262
1263
1264     for(Int_t nh = 0; nh < 9; nh++) {
1265
1266       if( eta > etaHigh[nh]/10.0 || eta < etaLow[nh]/10.0 )
1267         continue;
1268
1269
1270       if(fAnalysisMC){
1271         hMcOut[0][nh]->Fill(esdTrack->Pt());
1272         hMcOut[pidCode][nh]->Fill(esdTrack->Pt());
1273       }
1274
1275       histAllCh[nh]->Fill(momentum, dedx);
1276       if(beta>1){
1277         histPiTof[nh]->Fill(momentum, dedx);
1278         histpPiTof[nh]->Fill(momentum);
1279       }
1280
1281       if( momentum <= 0.6 && momentum >= 0.4  ){//only p:0.4-0.6 GeV, pion MIP
1282         //Fill  pion MIP, before calibration
1283         if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
1284           hMIPVsPhi[nh]->Fill(phi,dedx);
1285           pMIPVsPhi[nh]->Fill(phi,dedx);
1286
1287
1288           hMIPVsNch[nh]->Fill(multTPC,dedx);
1289           pMIPVsNch[nh]->Fill(multTPC,dedx);
1290
1291         }
1292
1293         //Fill electrons, before calibration
1294         if( dedx > DeDxMIPMax+1 && dedx < 95 ){
1295           if(TMath::Abs(beta-1)<0.1){
1296             hPlateauVsPhi[nh]->Fill(phi,dedx);
1297             pPlateauVsPhi[nh]->Fill(phi,dedx);
1298           }
1299         }
1300
1301       }
1302
1303
1304     }//end loop over eta intervals
1305
1306
1307
1308
1309
1310   }//end of track loop
1311
1312
1313
1314
1315 }
1316 //__________________________________________________________________
1317 void AliAnalysisTaskQAHighPtDeDx::ProduceArrayTrksAOD( AliAODEvent *AODevent ){
1318
1319
1320   Int_t nAODTracks = AODevent->GetNumberOfTracks();
1321   Int_t multTPC = 0;
1322
1323   //get multiplicity tpc only track cuts
1324   for(Int_t iT = 0; iT < nAODTracks; iT++) {
1325
1326     AliAODTrack* aodTrack = AODevent->GetTrack(iT);
1327
1328     if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
1329       continue;
1330
1331
1332     if (fTrackFilterTPC) {
1333       // TPC only cuts is bit 1
1334       if(!aodTrack->TestFilterBit(1))
1335         continue;
1336     }
1337
1338     multTPC++;
1339
1340   }
1341
1342
1343   for(Int_t iT = 0; iT < nAODTracks; iT++) {
1344
1345     AliAODTrack* aodTrack = AODevent->GetTrack(iT);
1346
1347     if (fTrackFilterGolden) {
1348       // "Global track RAA analysis QM2011 + Chi2ITS<36"; bit 1024
1349       if(!aodTrack->TestFilterBit(1024))
1350         continue;
1351     }
1352
1353
1354     if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
1355       continue;
1356
1357
1358     Double_t eta  = aodTrack->Eta();
1359     Double_t phi  = aodTrack->Phi();
1360     Double_t momentum = aodTrack->P();
1361
1362
1363     if(!PhiCut(aodTrack->Pt(), phi, aodTrack->Charge(), magf, cutLow, cutHigh))
1364       continue;
1365
1366
1367
1368     AliAODPid* aodPid = aodTrack->GetDetPid();
1369     Short_t ncl     = -10;
1370     Float_t dedx    = -10;
1371
1372     //TOF
1373     Float_t beta = -99;
1374
1375
1376     if(aodPid) {
1377       ncl     = aodPid->GetTPCsignalN();
1378       dedx    = aodPid->GetTPCsignal();
1379       //TOF
1380       Bool_t IsTOFout=kFALSE;
1381       Float_t lengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF
1382       Float_t timeTOF=-999;
1383
1384       if ((aodTrack->GetStatus()&AliESDtrack::kTOFout)==0)
1385         IsTOFout=kTRUE;
1386
1387       lengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF
1388
1389       timeTOF=aodPid->GetTOFsignal();
1390
1391       Double_t inttime[5]={0,0,0,0,0};
1392       aodPid->GetIntegratedTimes(inttime);// Returns the array with integrated times for each particle hypothesis
1393
1394
1395       if ( !IsTOFout ){
1396         if ( ( lengthtrack != 0 ) && ( timeTOF != 0) )
1397           beta = inttime[0] / timeTOF;
1398       }
1399
1400     }
1401
1402
1403     if(ncl<70)
1404       continue;
1405
1406
1407     Short_t pidCode     = 0;
1408
1409     if(fAnalysisMC) {
1410
1411       const Int_t label = TMath::Abs(aodTrack->GetLabel());
1412       AliAODMCParticle* mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(label));
1413
1414       if (mcTrack){
1415
1416         Int_t pdgCode = mcTrack->GetPdgCode();
1417         pidCode = GetPidCode(pdgCode);
1418
1419       }
1420
1421     }
1422
1423
1424
1425
1426     //if(momentum>0.6||momentum<0.4)//only p:0.4-0.6 GeV
1427         //continue;
1428
1429     //etaLow
1430     //etaHigh
1431     if( momentum <= 0.6 && momentum >= 0.4  ){//only p:0.4-0.6 GeV, pion MIP
1432       if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
1433         if(momentum<0.6&&momentum>0.4){
1434           hMIPVsEta->Fill(eta,dedx);
1435           pMIPVsEta->Fill(eta,dedx);
1436         }
1437       }
1438       if( dedx > DeDxMIPMax+1 && dedx < 95 ){
1439         if(TMath::Abs(beta-1)<0.1){
1440           hPlateauVsEta->Fill(eta,dedx);
1441           pPlateauVsEta->Fill(eta,dedx);
1442         }
1443       }
1444     }
1445
1446
1447     for(Int_t nh = 0; nh < 9; nh++) {
1448
1449       if( eta > etaHigh[nh]/10.0 || eta < etaLow[nh]/10.0 )
1450         continue;
1451
1452
1453       if(fAnalysisMC){
1454         hMcOut[0][nh]->Fill(aodTrack->Pt());
1455         hMcOut[pidCode][nh]->Fill(aodTrack->Pt());
1456       }
1457
1458       histAllCh[nh]->Fill(momentum, dedx);
1459       if(beta>1){
1460         histPiTof[nh]->Fill(momentum, dedx);
1461         histpPiTof[nh]->Fill(momentum);
1462       }
1463
1464       if( momentum <= 0.6 && momentum >= 0.4  ){//only p:0.4-0.6 GeV, pion MIP
1465         //Fill  pion MIP, before calibration
1466         if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
1467           hMIPVsPhi[nh]->Fill(phi,dedx);
1468           pMIPVsPhi[nh]->Fill(phi,dedx);
1469
1470
1471           hMIPVsNch[nh]->Fill(multTPC,dedx);
1472           pMIPVsNch[nh]->Fill(multTPC,dedx);
1473
1474         }
1475
1476         //Fill electrons, before calibration
1477         if( dedx > DeDxMIPMax+1 && dedx < 95 ){
1478           if(TMath::Abs(beta-1)<0.1){
1479             hPlateauVsPhi[nh]->Fill(phi,dedx);
1480             pPlateauVsPhi[nh]->Fill(phi,dedx);
1481           }
1482         }
1483
1484       }
1485
1486     }//end loop over eta intervals
1487
1488
1489
1490
1491
1492   }//end of track loop
1493
1494
1495
1496
1497
1498
1499
1500 }
1501 //__________________________________________________
1502 void AliAnalysisTaskQAHighPtDeDx::ProduceArrayV0ESD( AliESDEvent *ESDevent ){
1503   Int_t nv0s = ESDevent->GetNumberOfV0s();
1504   /*
1505   if(nv0s<1){
1506     return;
1507     }*/
1508
1509   const AliESDVertex *myBestPrimaryVertex = ESDevent->GetPrimaryVertex();
1510   if (!myBestPrimaryVertex) return;
1511   if (!(myBestPrimaryVertex->GetStatus())) return;
1512   Double_t  lPrimaryVtxPosition[3];
1513   myBestPrimaryVertex->GetXYZ(lPrimaryVtxPosition);
1514
1515   Double_t  lPrimaryVtxCov[6];
1516   myBestPrimaryVertex->GetCovMatrix(lPrimaryVtxCov);
1517   Double_t  lPrimaryVtxChi2 = myBestPrimaryVertex->GetChi2toNDF();
1518
1519   AliAODVertex* myPrimaryVertex = new AliAODVertex(lPrimaryVtxPosition, lPrimaryVtxCov, lPrimaryVtxChi2, NULL, -1, AliAODVertex::kPrimary);
1520
1521
1522
1523   //
1524   // LOOP OVER V0s, K0s, L, AL
1525   //
1526
1527
1528   for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
1529
1530     // This is the begining of the V0 loop
1531     AliESDv0 *esdV0 = ESDevent->GetV0(iV0);
1532     if (!esdV0) continue;
1533
1534     //check onfly status
1535     if(esdV0->GetOnFlyStatus()!=0)
1536       continue;
1537
1538
1539     // AliESDTrack (V0 Daughters)
1540     UInt_t lKeyPos = (UInt_t)TMath::Abs(esdV0->GetPindex());
1541     UInt_t lKeyNeg = (UInt_t)TMath::Abs(esdV0->GetNindex());
1542
1543     AliESDtrack *pTrack = ESDevent->GetTrack(lKeyPos);
1544     AliESDtrack *nTrack = ESDevent->GetTrack(lKeyNeg);
1545     if (!pTrack || !nTrack) {
1546       Printf("ERROR: Could not retreive one of the daughter track");
1547       continue;
1548     }
1549
1550     // Remove like-sign
1551     if (pTrack->GetSign() == nTrack->GetSign()) {
1552       //cout<< "like sign, continue"<< endl;
1553       continue;
1554     }
1555
1556     // Eta cut on decay products
1557     if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
1558       continue;
1559
1560
1561     // Check if switch does anything!
1562     Bool_t isSwitched = kFALSE;
1563     if (pTrack->GetSign() < 0) { // switch
1564
1565       isSwitched = kTRUE;
1566       AliESDtrack* helpTrack = nTrack;
1567       nTrack = pTrack;
1568       pTrack = helpTrack;
1569     }
1570
1571     //Double_t  lV0Position[3];
1572     //esdV0->GetXYZ(lV0Position[0], lV0Position[1], lV0Position[2]);
1573
1574
1575     AliKFVertex primaryVtxKF( *myPrimaryVertex );
1576     AliKFParticle::SetField(ESDevent->GetMagneticField());
1577
1578     // Also implement switch here!!!!!!
1579     AliKFParticle* negEKF  = 0; // e-
1580     AliKFParticle* posEKF  = 0; // e+
1581     AliKFParticle* negPiKF = 0; // pi -
1582     AliKFParticle* posPiKF = 0; // pi +
1583     AliKFParticle* posPKF  = 0; // p
1584     AliKFParticle* negAPKF = 0; // p-bar
1585
1586     if(!isSwitched) {
1587       negEKF  = new AliKFParticle( *(esdV0->GetParamN()) , 11);
1588       posEKF  = new AliKFParticle( *(esdV0->GetParamP()) ,-11);
1589       negPiKF = new AliKFParticle( *(esdV0->GetParamN()) ,-211);
1590       posPiKF = new AliKFParticle( *(esdV0->GetParamP()) , 211);
1591       posPKF  = new AliKFParticle( *(esdV0->GetParamP()) , 2212);
1592       negAPKF = new AliKFParticle( *(esdV0->GetParamN()) ,-2212);
1593     } else { // switch + and -
1594       negEKF  = new AliKFParticle( *(esdV0->GetParamP()) , 11);
1595       posEKF  = new AliKFParticle( *(esdV0->GetParamN()) ,-11);
1596       negPiKF = new AliKFParticle( *(esdV0->GetParamP()) ,-211);
1597       posPiKF = new AliKFParticle( *(esdV0->GetParamN()) , 211);
1598       posPKF  = new AliKFParticle( *(esdV0->GetParamN()) , 2212);
1599       negAPKF = new AliKFParticle( *(esdV0->GetParamP()) ,-2212);
1600     }
1601
1602     AliKFParticle v0GKF;  // Gamma e.g. from pi0
1603     v0GKF+=(*negEKF);
1604     v0GKF+=(*posEKF);
1605     v0GKF.SetProductionVertex(primaryVtxKF);
1606
1607     AliKFParticle v0K0sKF; // K0 short
1608     v0K0sKF+=(*negPiKF);
1609     v0K0sKF+=(*posPiKF);
1610     v0K0sKF.SetProductionVertex(primaryVtxKF);
1611
1612     AliKFParticle v0LambdaKF; // Lambda
1613     v0LambdaKF+=(*negPiKF);
1614     v0LambdaKF+=(*posPKF);
1615     v0LambdaKF.SetProductionVertex(primaryVtxKF);
1616
1617     AliKFParticle v0AntiLambdaKF; // Lambda-bar
1618     v0AntiLambdaKF+=(*posPiKF);
1619     v0AntiLambdaKF+=(*negAPKF);
1620     v0AntiLambdaKF.SetProductionVertex(primaryVtxKF);
1621
1622     Double_t dmassG     = v0GKF.GetMass();
1623     Double_t dmassK     = v0K0sKF.GetMass()-0.498;
1624     Double_t dmassL     = v0LambdaKF.GetMass()-1.116;
1625     Double_t dmassAL    = v0AntiLambdaKF.GetMass()-1.116;
1626
1627
1628     for( Int_t case_v0 = 0; case_v0 < 2; ++case_v0 ){
1629
1630
1631       switch(case_v0){
1632       case 0:{
1633
1634         Bool_t fillPos = kFALSE;
1635         Bool_t fillNeg = kFALSE;
1636
1637         if(dmassG < 0.1)
1638           continue;
1639
1640         if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01){
1641           continue;
1642         }
1643
1644         if(dmassL<0.01){
1645           fillPos = kTRUE;
1646         }
1647         if(dmassAL<0.01) {
1648           if(fillPos)
1649             continue;
1650           fillNeg = kTRUE;
1651         }
1652
1653         if(dmassK<0.01) {
1654           if(fillPos||fillNeg)
1655             continue;
1656           fillPos = kTRUE;
1657           fillNeg = kTRUE;
1658         }
1659
1660
1661         for(Int_t j = 0; j < 2; j++) {
1662
1663           AliESDtrack* track = 0;
1664
1665           if(j==0) {
1666
1667             if(fillNeg)
1668               track = nTrack;
1669             else
1670               continue;
1671           } else {
1672
1673             if(fillPos)
1674               track = pTrack;
1675             else
1676               continue;
1677           }
1678
1679           if(track->GetTPCsignalN()<=70)continue;
1680           Double_t phi     = track->Phi();
1681
1682           if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
1683             continue;
1684
1685
1686           //if(!PhiCut(pt, phi, charge, magf, cutLow, cutHigh, hPhi))
1687           //    continue;
1688
1689           Double_t eta  = track->Eta();
1690           Double_t momentum = track->Pt();
1691           Double_t dedx = track->GetTPCsignal();
1692
1693           if(fillPos&&fillNeg){
1694
1695
1696             if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
1697               if(momentum<0.6&&momentum>0.4){
1698                 hMIPVsEtaV0s->Fill(eta,dedx);
1699                 pMIPVsEtaV0s->Fill(eta,dedx);
1700               }
1701             }
1702
1703
1704           }
1705
1706           for(Int_t nh = 0; nh < nHists; nh++) {
1707
1708
1709
1710             if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
1711               continue;
1712
1713             if(fillPos&&fillNeg){
1714
1715               histPiV0[nh]->Fill(momentum, dedx);
1716               histpPiV0[nh]->Fill(momentum);
1717
1718             }
1719             else{
1720
1721               histPV0[nh]->Fill(momentum, dedx);
1722               histpPV0[nh]->Fill(momentum);
1723
1724             }
1725
1726           }
1727
1728         }//end loop over two tracks
1729
1730       };
1731         break;
1732
1733       case 1:{//gammas
1734
1735         Bool_t fillPos = kFALSE;
1736         Bool_t fillNeg = kFALSE;
1737
1738
1739
1740         if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01) {
1741           if(dmassG<0.01 && dmassG>0.0001) {
1742
1743
1744             if( TMath::Abs(nTrack->GetTPCsignal() - 85.0) < 5)
1745               fillPos = kTRUE;
1746             if( TMath::Abs(pTrack->GetTPCsignal() - 85.0) < 5)
1747               fillNeg = kTRUE;
1748
1749           } else {
1750             continue;
1751           }
1752         }
1753
1754         //cout<<"fillPos="<<fillPos<<endl;
1755
1756         if(fillPos == kTRUE && fillNeg == kTRUE)
1757           continue;
1758
1759
1760         AliESDtrack* track = 0;
1761         if(fillNeg)
1762           track = nTrack;
1763         else if(fillPos)
1764           track = pTrack;
1765         else
1766           continue;
1767
1768         Double_t dedx  = track->GetTPCsignal();
1769         Double_t eta  = track->Eta();
1770         Double_t phi  = track->Phi();
1771         Double_t momentum = track->P();
1772
1773         if(track->GetTPCsignalN()<=70)continue;
1774
1775         if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
1776           continue;
1777
1778         for(Int_t nh = 0; nh < nHists; nh++) {
1779
1780           if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
1781             continue;
1782
1783           histEV0[nh]->Fill(momentum, dedx);
1784
1785         }
1786
1787       };
1788         break;
1789
1790
1791       }//end switch
1792
1793     }//end loop over case V0
1794
1795
1796     // clean up loop over v0
1797
1798     delete negPiKF;
1799     delete posPiKF;
1800     delete posPKF;
1801     delete negAPKF;
1802
1803
1804
1805   }
1806
1807
1808   delete myPrimaryVertex;
1809
1810
1811 }
1812 //__________________________________________________________________________
1813 void AliAnalysisTaskQAHighPtDeDx::ProduceArrayV0AOD( AliAODEvent *AODevent ){
1814   Int_t nv0s = AODevent->GetNumberOfV0s();
1815   /*
1816   if(nv0s<1){
1817     return;
1818     }*/
1819
1820   AliAODVertex *myBestPrimaryVertex = AODevent->GetPrimaryVertex();
1821   if (!myBestPrimaryVertex) return;
1822
1823
1824
1825   // ################################
1826   // #### BEGINNING OF V0 CODE ######
1827   // ################################
1828   // This is the begining of the V0 loop
1829   for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
1830     AliAODv0 *aodV0 = AODevent->GetV0(iV0);
1831     if (!aodV0) continue;
1832
1833
1834     //check onfly status
1835     if(aodV0->GetOnFlyStatus()!=0)
1836       continue;
1837
1838     // AliAODTrack (V0 Daughters)
1839     AliAODVertex* vertex = aodV0->GetSecondaryVtx();
1840     if (!vertex) {
1841       Printf("ERROR: Could not retrieve vertex");
1842       continue;
1843     }
1844
1845     AliAODTrack *pTrack = (AliAODTrack*)vertex->GetDaughter(0);
1846     AliAODTrack *nTrack = (AliAODTrack*)vertex->GetDaughter(1);
1847     if (!pTrack || !nTrack) {
1848       Printf("ERROR: Could not retrieve one of the daughter track");
1849       continue;
1850     }
1851
1852     // Remove like-sign
1853     if (pTrack->Charge() == nTrack->Charge()) {
1854       //cout<< "like sign, continue"<< endl;
1855       continue;
1856     }
1857
1858     // Make sure charge ordering is ok
1859     if (pTrack->Charge() < 0) {
1860       AliAODTrack* helpTrack = pTrack;
1861       pTrack = nTrack;
1862       nTrack = helpTrack;
1863     }
1864
1865     // Eta cut on decay products
1866     if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
1867       continue;
1868
1869
1870     Double_t dmassG  = aodV0->InvMass2Prongs(0,1,11,11);
1871     Double_t dmassK  = aodV0->MassK0Short()-0.498;
1872     Double_t dmassL  = aodV0->MassLambda()-1.116;
1873     Double_t dmassAL = aodV0->MassAntiLambda()-1.116;
1874
1875     for( Int_t case_v0 = 0; case_v0 < 2; ++case_v0 ){
1876
1877
1878       switch(case_v0){
1879       case 0:{
1880         Bool_t fillPos = kFALSE;
1881         Bool_t fillNeg = kFALSE;
1882
1883
1884         if(dmassG < 0.1)
1885           continue;
1886
1887         if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01){
1888           continue;
1889         }
1890
1891         if(dmassL<0.01){
1892           fillPos = kTRUE;
1893         }
1894         if(dmassAL<0.01) {
1895           if(fillPos)
1896             continue;
1897           fillNeg = kTRUE;
1898         }
1899
1900         if(dmassK<0.01) {
1901           if(fillPos||fillNeg)
1902             continue;
1903           fillPos = kTRUE;
1904           fillNeg = kTRUE;
1905         }
1906
1907
1908
1909         for(Int_t j = 0; j < 2; j++) {
1910
1911           AliAODTrack* track = 0;
1912
1913           if(j==0) {
1914
1915             if(fillNeg)
1916               track = nTrack;
1917             else
1918               continue;
1919           } else {
1920
1921             if(fillPos)
1922               track = pTrack;
1923             else
1924               continue;
1925           }
1926
1927           if(track->GetTPCsignalN()<=70)continue;
1928
1929           Double_t phi     = track->Phi();
1930
1931           if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
1932             continue;
1933
1934
1935           //if(!PhiCut(pt, phi, charge, magf, cutLow, cutHigh, hPhi))
1936           //    continue;
1937
1938           Double_t eta  = track->Eta();
1939           Double_t momentum = track->Pt();
1940           Double_t dedx = track->GetTPCsignal();
1941
1942           if(fillPos&&fillNeg){
1943
1944
1945             if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
1946               if(momentum<0.6&&momentum>0.4){
1947                 hMIPVsEtaV0s->Fill(eta,dedx);
1948                 pMIPVsEtaV0s->Fill(eta,dedx);
1949               }
1950             }
1951
1952
1953           }
1954
1955           for(Int_t nh = 0; nh < nHists; nh++) {
1956
1957
1958
1959             if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
1960               continue;
1961
1962             if(fillPos&&fillNeg){
1963
1964               histPiV0[nh]->Fill(momentum, dedx);
1965               histpPiV0[nh]->Fill(momentum);
1966
1967             }
1968             else{
1969
1970               histPV0[nh]->Fill(momentum, dedx);
1971               histpPV0[nh]->Fill(momentum);
1972
1973             }
1974
1975           }
1976
1977
1978         }//end loop over two tracks
1979       };
1980         break;
1981
1982       case 1:{//gammas
1983
1984         Bool_t fillPos = kFALSE;
1985         Bool_t fillNeg = kFALSE;
1986
1987         if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01) {
1988           if(dmassG<0.01 && dmassG>0.0001) {
1989
1990             if( TMath::Abs(nTrack->GetTPCsignal() - 85.0) < 5)
1991               fillPos = kTRUE;
1992             if( TMath::Abs(pTrack->GetTPCsignal() - 85.0) < 5)
1993               fillNeg = kTRUE;
1994
1995           } else {
1996             continue;
1997           }
1998         }
1999
2000
2001         if(fillPos == kTRUE && fillNeg == kTRUE)
2002           continue;
2003
2004
2005         AliAODTrack* track = 0;
2006         if(fillNeg)
2007           track = nTrack;
2008         else if(fillPos)
2009           track = pTrack;
2010         else
2011           continue;
2012
2013         Double_t dedx  = track->GetTPCsignal();
2014         Double_t eta  = track->Eta();
2015         Double_t phi  = track->Phi();
2016         Double_t momentum = track->P();
2017
2018         if(track->GetTPCsignalN()<=70)continue;
2019
2020         if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
2021           continue;
2022
2023         for(Int_t nh = 0; nh < nHists; nh++) {
2024
2025           if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
2026             continue;
2027
2028           histEV0[nh]->Fill(momentum, dedx);
2029
2030         }
2031
2032       };
2033         break;
2034
2035
2036       }//end switch
2037     }//end loop over V0s cases
2038
2039   }//end loop over v0's
2040
2041
2042
2043
2044 }
2045 Bool_t AliAnalysisTaskQAHighPtDeDx::PhiCut(Double_t pt, Double_t phi, Double_t q, Float_t   mag, TF1* phiCutLow, TF1* phiCutHigh)
2046 {
2047   if(pt < 2.0)
2048     return kTRUE;
2049
2050   //Double_t phi = track->Phi();
2051   if(mag < 0)    // for negatve polarity field
2052     phi = TMath::TwoPi() - phi;
2053   if(q < 0) // for negatve charge
2054     phi = TMath::TwoPi()-phi;
2055
2056   phi += TMath::Pi()/18.0; // to center gap in the middle
2057   phi = fmod(phi, TMath::Pi()/9.0);
2058
2059   if(phi<phiCutHigh->Eval(pt)
2060      && phi>phiCutLow->Eval(pt))
2061   return kFALSE; // reject track
2062
2063   hPhi->Fill(pt, phi);
2064
2065   return kTRUE;
2066 }