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