coverity fix and some clean up
[u/mrichter/AliRoot.git] / 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 = dynamic_cast<AliAODTrack*>(AODevent->GetTrack(iT));
1327     if(!aodTrack) AliFatal("Not a standard AOD");
1328     
1329     if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
1330       continue;
1331     
1332
1333     if (fTrackFilterTPC) {
1334       // TPC only cuts is bit 1
1335       if(!aodTrack->TestFilterBit(1))
1336         continue;
1337     }
1338    
1339     multTPC++;
1340     
1341   }
1342
1343
1344   for(Int_t iT = 0; iT < nAODTracks; iT++) {
1345     
1346     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*>(AODevent->GetTrack(iT));
1347     if(!aodTrack) AliFatal("Not a standard AOD");
1348     
1349     if (fTrackFilterGolden) {     
1350       // "Global track RAA analysis QM2011 + Chi2ITS<36"; bit 1024
1351       if(!aodTrack->TestFilterBit(1024))
1352         continue;
1353     }
1354     
1355     
1356     if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
1357       continue;
1358     
1359     
1360     Double_t eta  = aodTrack->Eta();
1361     Double_t phi  = aodTrack->Phi();
1362     Double_t momentum = aodTrack->P();
1363     
1364  
1365     if(!PhiCut(aodTrack->Pt(), phi, aodTrack->Charge(), magf, cutLow, cutHigh))
1366       continue;
1367     
1368
1369     
1370     AliAODPid* aodPid = aodTrack->GetDetPid();
1371     Short_t ncl     = -10;
1372     Float_t dedx    = -10;
1373
1374     //TOF    
1375     Float_t beta = -99;
1376
1377
1378     if(aodPid) {
1379       ncl     = aodPid->GetTPCsignalN();
1380       dedx    = aodPid->GetTPCsignal();
1381       //TOF
1382       Bool_t IsTOFout=kFALSE;
1383       Float_t lengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF 
1384       Float_t timeTOF=-999;
1385
1386       if ((aodTrack->GetStatus()&AliESDtrack::kTOFout)==0)
1387         IsTOFout=kTRUE;
1388       
1389       lengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF 
1390       
1391       timeTOF=aodPid->GetTOFsignal();
1392       
1393       Double_t inttime[5]={0,0,0,0,0}; 
1394       aodPid->GetIntegratedTimes(inttime);// Returns the array with integrated times for each particle hypothesis
1395
1396       
1397       if ( !IsTOFout ){
1398         if ( ( lengthtrack != 0 ) && ( timeTOF != 0) )
1399           beta = inttime[0] / timeTOF;
1400       }
1401       
1402     }
1403
1404
1405     if(ncl<70)
1406       continue;
1407
1408
1409     Short_t pidCode     = 0; 
1410     
1411     if(fAnalysisMC) {
1412       
1413       const Int_t label = TMath::Abs(aodTrack->GetLabel());
1414       AliAODMCParticle* mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(label));
1415       
1416       if (mcTrack){
1417         
1418         Int_t pdgCode = mcTrack->GetPdgCode();
1419         pidCode = GetPidCode(pdgCode);
1420         
1421       }
1422       
1423     }
1424
1425
1426
1427
1428     //if(momentum>0.6||momentum<0.4)//only p:0.4-0.6 GeV
1429         //continue;
1430     
1431     //etaLow
1432     //etaHigh
1433     if( momentum <= 0.6 && momentum >= 0.4  ){//only p:0.4-0.6 GeV, pion MIP
1434       if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){       
1435         if(momentum<0.6&&momentum>0.4){
1436           hMIPVsEta->Fill(eta,dedx);
1437           pMIPVsEta->Fill(eta,dedx);
1438         }
1439       }
1440       if( dedx > DeDxMIPMax+1 && dedx < 95 ){
1441         if(TMath::Abs(beta-1)<0.1){
1442           hPlateauVsEta->Fill(eta,dedx);
1443           pPlateauVsEta->Fill(eta,dedx); 
1444         }
1445       }
1446     }
1447     
1448     
1449     for(Int_t nh = 0; nh < 9; nh++) {
1450       
1451       if( eta > etaHigh[nh]/10.0 || eta < etaLow[nh]/10.0 )
1452         continue;
1453       
1454       
1455       if(fAnalysisMC){
1456         hMcOut[0][nh]->Fill(aodTrack->Pt());
1457         hMcOut[pidCode][nh]->Fill(aodTrack->Pt());
1458       }
1459       
1460       histAllCh[nh]->Fill(momentum, dedx);
1461       if(beta>1){
1462         histPiTof[nh]->Fill(momentum, dedx);
1463         histpPiTof[nh]->Fill(momentum);
1464       }
1465
1466       if( momentum <= 0.6 && momentum >= 0.4  ){//only p:0.4-0.6 GeV, pion MIP
1467         //Fill  pion MIP, before calibration
1468         if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){     
1469           hMIPVsPhi[nh]->Fill(phi,dedx);
1470           pMIPVsPhi[nh]->Fill(phi,dedx);
1471           
1472           
1473           hMIPVsNch[nh]->Fill(multTPC,dedx);
1474           pMIPVsNch[nh]->Fill(multTPC,dedx);
1475           
1476         }
1477         
1478         //Fill electrons, before calibration
1479         if( dedx > DeDxMIPMax+1 && dedx < 95 ){
1480           if(TMath::Abs(beta-1)<0.1){
1481             hPlateauVsPhi[nh]->Fill(phi,dedx);
1482             pPlateauVsPhi[nh]->Fill(phi,dedx);
1483           }
1484         }
1485
1486       }
1487
1488     }//end loop over eta intervals
1489     
1490     
1491
1492     
1493     
1494   }//end of track loop
1495   
1496  
1497   
1498   
1499   
1500  
1501   
1502 }
1503 //__________________________________________________
1504 void AliAnalysisTaskQAHighPtDeDx::ProduceArrayV0ESD( AliESDEvent *ESDevent ){
1505   Int_t nv0s = ESDevent->GetNumberOfV0s();
1506   /*
1507   if(nv0s<1){
1508     return;
1509     }*/
1510   
1511   const AliESDVertex *myBestPrimaryVertex = ESDevent->GetPrimaryVertex();
1512   if (!myBestPrimaryVertex) return;
1513   if (!(myBestPrimaryVertex->GetStatus())) return;
1514   Double_t  lPrimaryVtxPosition[3];
1515   myBestPrimaryVertex->GetXYZ(lPrimaryVtxPosition);
1516   
1517   Double_t  lPrimaryVtxCov[6];
1518   myBestPrimaryVertex->GetCovMatrix(lPrimaryVtxCov);
1519   Double_t  lPrimaryVtxChi2 = myBestPrimaryVertex->GetChi2toNDF();
1520   
1521   AliAODVertex* myPrimaryVertex = new AliAODVertex(lPrimaryVtxPosition, lPrimaryVtxCov, lPrimaryVtxChi2, NULL, -1, AliAODVertex::kPrimary);
1522
1523
1524     
1525   //
1526   // LOOP OVER V0s, K0s, L, AL
1527   //
1528   
1529   
1530   for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
1531     
1532     // This is the begining of the V0 loop  
1533     AliESDv0 *esdV0 = ESDevent->GetV0(iV0);
1534     if (!esdV0) continue;
1535     
1536     //check onfly status
1537     if(esdV0->GetOnFlyStatus()!=0)
1538       continue;
1539
1540
1541     // AliESDTrack (V0 Daughters)
1542     UInt_t lKeyPos = (UInt_t)TMath::Abs(esdV0->GetPindex());
1543     UInt_t lKeyNeg = (UInt_t)TMath::Abs(esdV0->GetNindex());
1544     
1545     AliESDtrack *pTrack = ESDevent->GetTrack(lKeyPos);
1546     AliESDtrack *nTrack = ESDevent->GetTrack(lKeyNeg);
1547     if (!pTrack || !nTrack) {
1548       Printf("ERROR: Could not retreive one of the daughter track");
1549       continue;
1550     }
1551     
1552     // Remove like-sign
1553     if (pTrack->GetSign() == nTrack->GetSign()) {
1554       //cout<< "like sign, continue"<< endl;
1555       continue;
1556     } 
1557     
1558     // Eta cut on decay products
1559     if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
1560       continue;
1561     
1562     
1563     // Check if switch does anything!
1564     Bool_t isSwitched = kFALSE;
1565     if (pTrack->GetSign() < 0) { // switch
1566       
1567       isSwitched = kTRUE;
1568       AliESDtrack* helpTrack = nTrack;
1569       nTrack = pTrack;
1570       pTrack = helpTrack;
1571     }   
1572     
1573     //Double_t  lV0Position[3];
1574     //esdV0->GetXYZ(lV0Position[0], lV0Position[1], lV0Position[2]);
1575     
1576
1577     AliKFVertex primaryVtxKF( *myPrimaryVertex );
1578     AliKFParticle::SetField(ESDevent->GetMagneticField());
1579     
1580     // Also implement switch here!!!!!!
1581     AliKFParticle* negEKF  = 0; // e-
1582     AliKFParticle* posEKF  = 0; // e+
1583     AliKFParticle* negPiKF = 0; // pi -
1584     AliKFParticle* posPiKF = 0; // pi +
1585     AliKFParticle* posPKF  = 0; // p
1586     AliKFParticle* negAPKF = 0; // p-bar
1587     
1588     if(!isSwitched) {
1589       negEKF  = new AliKFParticle( *(esdV0->GetParamN()) , 11);
1590       posEKF  = new AliKFParticle( *(esdV0->GetParamP()) ,-11);
1591       negPiKF = new AliKFParticle( *(esdV0->GetParamN()) ,-211);
1592       posPiKF = new AliKFParticle( *(esdV0->GetParamP()) , 211);
1593       posPKF  = new AliKFParticle( *(esdV0->GetParamP()) , 2212);
1594       negAPKF = new AliKFParticle( *(esdV0->GetParamN()) ,-2212);
1595     } else { // switch + and - 
1596       negEKF  = new AliKFParticle( *(esdV0->GetParamP()) , 11);
1597       posEKF  = new AliKFParticle( *(esdV0->GetParamN()) ,-11);
1598       negPiKF = new AliKFParticle( *(esdV0->GetParamP()) ,-211);
1599       posPiKF = new AliKFParticle( *(esdV0->GetParamN()) , 211);
1600       posPKF  = new AliKFParticle( *(esdV0->GetParamN()) , 2212);
1601       negAPKF = new AliKFParticle( *(esdV0->GetParamP()) ,-2212);
1602     }
1603     
1604     AliKFParticle v0GKF;  // Gamma e.g. from pi0
1605     v0GKF+=(*negEKF);
1606     v0GKF+=(*posEKF);
1607     v0GKF.SetProductionVertex(primaryVtxKF);
1608     
1609     AliKFParticle v0K0sKF; // K0 short
1610     v0K0sKF+=(*negPiKF);
1611     v0K0sKF+=(*posPiKF);
1612     v0K0sKF.SetProductionVertex(primaryVtxKF);
1613     
1614     AliKFParticle v0LambdaKF; // Lambda
1615     v0LambdaKF+=(*negPiKF);
1616     v0LambdaKF+=(*posPKF);      
1617     v0LambdaKF.SetProductionVertex(primaryVtxKF);
1618     
1619     AliKFParticle v0AntiLambdaKF; // Lambda-bar
1620     v0AntiLambdaKF+=(*posPiKF);
1621     v0AntiLambdaKF+=(*negAPKF);
1622     v0AntiLambdaKF.SetProductionVertex(primaryVtxKF);
1623     
1624     Double_t dmassG     = v0GKF.GetMass();
1625     Double_t dmassK     = v0K0sKF.GetMass()-0.498;
1626     Double_t dmassL     = v0LambdaKF.GetMass()-1.116;
1627     Double_t dmassAL    = v0AntiLambdaKF.GetMass()-1.116;
1628
1629
1630     for( Int_t case_v0 = 0; case_v0 < 2; ++case_v0 ){
1631
1632
1633       switch(case_v0){
1634       case 0:{    
1635         
1636         Bool_t fillPos = kFALSE;
1637         Bool_t fillNeg = kFALSE;
1638         
1639         if(dmassG < 0.1)
1640           continue;
1641         
1642         if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01){
1643           continue;
1644         }
1645         
1646         if(dmassL<0.01){
1647           fillPos = kTRUE;
1648         }
1649         if(dmassAL<0.01) {
1650           if(fillPos)
1651             continue;
1652           fillNeg = kTRUE;
1653         }
1654         
1655         if(dmassK<0.01) {
1656           if(fillPos||fillNeg)
1657             continue;
1658           fillPos = kTRUE;
1659           fillNeg = kTRUE;
1660         }
1661         
1662         
1663         for(Int_t j = 0; j < 2; j++) {
1664           
1665           AliESDtrack* track = 0;
1666           
1667           if(j==0) {
1668             
1669             if(fillNeg)
1670               track = nTrack;
1671             else
1672               continue;
1673           } else {
1674             
1675             if(fillPos)
1676               track = pTrack;
1677             else
1678               continue;
1679           }
1680           
1681           if(track->GetTPCsignalN()<=70)continue;
1682           Double_t phi     = track->Phi();
1683           
1684           if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
1685             continue;
1686           
1687           
1688           //if(!PhiCut(pt, phi, charge, magf, cutLow, cutHigh, hPhi))
1689           //    continue;
1690           
1691           Double_t eta  = track->Eta();
1692           Double_t momentum = track->Pt();
1693           Double_t dedx = track->GetTPCsignal();
1694           
1695           if(fillPos&&fillNeg){
1696             
1697             
1698             if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){         
1699               if(momentum<0.6&&momentum>0.4){
1700                 hMIPVsEtaV0s->Fill(eta,dedx);
1701                 pMIPVsEtaV0s->Fill(eta,dedx);
1702               }
1703             }
1704             
1705             
1706           }
1707           
1708           for(Int_t nh = 0; nh < nHists; nh++) {
1709             
1710             
1711             
1712             if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
1713               continue;
1714             
1715             if(fillPos&&fillNeg){
1716               
1717               histPiV0[nh]->Fill(momentum, dedx);           
1718               histpPiV0[nh]->Fill(momentum);
1719
1720             }
1721             else{
1722               
1723               histPV0[nh]->Fill(momentum, dedx);
1724               histpPV0[nh]->Fill(momentum);
1725
1726             }
1727             
1728           }
1729                   
1730         }//end loop over two tracks
1731
1732       };
1733         break;
1734
1735       case 1:{//gammas    
1736         
1737         Bool_t fillPos = kFALSE;
1738         Bool_t fillNeg = kFALSE;
1739         
1740
1741
1742         if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01) {
1743           if(dmassG<0.01 && dmassG>0.0001) {
1744             
1745
1746             if( TMath::Abs(nTrack->GetTPCsignal() - 85.0) < 5)
1747               fillPos = kTRUE;
1748             if( TMath::Abs(pTrack->GetTPCsignal() - 85.0) < 5)
1749               fillNeg = kTRUE;
1750             
1751           } else {
1752             continue;
1753           }
1754         }
1755
1756         //cout<<"fillPos="<<fillPos<<endl;
1757
1758         if(fillPos == kTRUE && fillNeg == kTRUE)      
1759           continue;
1760         
1761         
1762         AliESDtrack* track = 0;
1763         if(fillNeg)
1764           track = nTrack;
1765         else if(fillPos)
1766           track = pTrack;
1767         else
1768           continue;
1769
1770         Double_t dedx  = track->GetTPCsignal();
1771         Double_t eta  = track->Eta();
1772         Double_t phi  = track->Phi();
1773         Double_t momentum = track->P();
1774
1775         if(track->GetTPCsignalN()<=70)continue;
1776
1777         if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
1778           continue;
1779         
1780         for(Int_t nh = 0; nh < nHists; nh++) {
1781       
1782           if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
1783             continue;
1784
1785           histEV0[nh]->Fill(momentum, dedx);
1786
1787         }
1788         
1789       };
1790         break;
1791         
1792         
1793       }//end switch
1794       
1795     }//end loop over case V0
1796
1797     
1798     // clean up loop over v0
1799     
1800     delete negPiKF;
1801     delete posPiKF;
1802     delete posPKF;
1803     delete negAPKF;
1804     
1805     
1806     
1807   }
1808   
1809   
1810   delete myPrimaryVertex;
1811   
1812   
1813 }
1814 //__________________________________________________________________________
1815 void AliAnalysisTaskQAHighPtDeDx::ProduceArrayV0AOD( AliAODEvent *AODevent ){
1816   Int_t nv0s = AODevent->GetNumberOfV0s();
1817   /*
1818   if(nv0s<1){
1819     return;
1820     }*/
1821   
1822   AliAODVertex *myBestPrimaryVertex = AODevent->GetPrimaryVertex();
1823   if (!myBestPrimaryVertex) return;
1824   
1825
1826      
1827   // ################################
1828   // #### BEGINNING OF V0 CODE ######
1829   // ################################
1830   // This is the begining of the V0 loop  
1831   for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
1832     AliAODv0 *aodV0 = AODevent->GetV0(iV0);
1833     if (!aodV0) continue;
1834     
1835     
1836     //check onfly status
1837     if(aodV0->GetOnFlyStatus()!=0)
1838       continue;
1839
1840     // AliAODTrack (V0 Daughters)
1841     AliAODVertex* vertex = aodV0->GetSecondaryVtx();
1842     if (!vertex) {
1843       Printf("ERROR: Could not retrieve vertex");
1844       continue;
1845     }
1846     
1847     AliAODTrack *pTrack = (AliAODTrack*)vertex->GetDaughter(0);
1848     AliAODTrack *nTrack = (AliAODTrack*)vertex->GetDaughter(1);
1849     if (!pTrack || !nTrack) {
1850       Printf("ERROR: Could not retrieve one of the daughter track");
1851       continue;
1852     }
1853     
1854     // Remove like-sign
1855     if (pTrack->Charge() == nTrack->Charge()) {
1856       //cout<< "like sign, continue"<< endl;
1857       continue;
1858     } 
1859     
1860     // Make sure charge ordering is ok
1861     if (pTrack->Charge() < 0) {
1862       AliAODTrack* helpTrack = pTrack;
1863       pTrack = nTrack;
1864       nTrack = helpTrack;
1865     } 
1866     
1867     // Eta cut on decay products
1868     if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
1869       continue;
1870     
1871     
1872     Double_t dmassG  = aodV0->InvMass2Prongs(0,1,11,11);
1873     Double_t dmassK  = aodV0->MassK0Short()-0.498;
1874     Double_t dmassL  = aodV0->MassLambda()-1.116;
1875     Double_t dmassAL = aodV0->MassAntiLambda()-1.116;
1876     
1877     for( Int_t case_v0 = 0; case_v0 < 2; ++case_v0 ){
1878       
1879       
1880       switch(case_v0){
1881       case 0:{   
1882         Bool_t fillPos = kFALSE;
1883         Bool_t fillNeg = kFALSE;
1884         
1885         
1886         if(dmassG < 0.1)
1887           continue;
1888         
1889         if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01){
1890           continue;
1891         }
1892         
1893         if(dmassL<0.01){
1894           fillPos = kTRUE;
1895         }
1896         if(dmassAL<0.01) {
1897           if(fillPos)
1898             continue;
1899           fillNeg = kTRUE;
1900         }
1901         
1902         if(dmassK<0.01) {
1903           if(fillPos||fillNeg)
1904             continue;
1905           fillPos = kTRUE;
1906           fillNeg = kTRUE;
1907         }
1908         
1909         
1910         
1911         for(Int_t j = 0; j < 2; j++) {
1912           
1913           AliAODTrack* track = 0;
1914           
1915           if(j==0) {
1916             
1917             if(fillNeg)
1918               track = nTrack;
1919             else
1920               continue;
1921           } else {
1922             
1923             if(fillPos)
1924               track = pTrack;
1925             else
1926               continue;
1927           }
1928           
1929           if(track->GetTPCsignalN()<=70)continue;
1930           
1931           Double_t phi     = track->Phi();
1932           
1933           if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
1934             continue;
1935           
1936           
1937           //if(!PhiCut(pt, phi, charge, magf, cutLow, cutHigh, hPhi))
1938           //    continue;
1939           
1940           Double_t eta  = track->Eta();
1941           Double_t momentum = track->Pt();
1942           Double_t dedx = track->GetTPCsignal();
1943           
1944           if(fillPos&&fillNeg){
1945             
1946             
1947             if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){         
1948               if(momentum<0.6&&momentum>0.4){
1949                 hMIPVsEtaV0s->Fill(eta,dedx);
1950                 pMIPVsEtaV0s->Fill(eta,dedx);
1951               }
1952             }
1953             
1954             
1955           }
1956           
1957           for(Int_t nh = 0; nh < nHists; nh++) {
1958             
1959             
1960             
1961             if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
1962               continue;
1963             
1964             if(fillPos&&fillNeg){
1965               
1966               histPiV0[nh]->Fill(momentum, dedx);           
1967               histpPiV0[nh]->Fill(momentum);        
1968
1969             }
1970             else{
1971               
1972               histPV0[nh]->Fill(momentum, dedx);
1973               histpPV0[nh]->Fill(momentum);
1974
1975             }
1976             
1977           }
1978           
1979
1980         }//end loop over two tracks
1981       };
1982         break;
1983         
1984       case 1:{//gammas    
1985         
1986         Bool_t fillPos = kFALSE;
1987         Bool_t fillNeg = kFALSE;
1988         
1989         if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01) {
1990           if(dmassG<0.01 && dmassG>0.0001) {
1991             
1992             if( TMath::Abs(nTrack->GetTPCsignal() - 85.0) < 5)
1993               fillPos = kTRUE;
1994             if( TMath::Abs(pTrack->GetTPCsignal() - 85.0) < 5)
1995               fillNeg = kTRUE;
1996             
1997           } else {
1998             continue;
1999           }
2000         }
2001
2002         
2003         if(fillPos == kTRUE && fillNeg == kTRUE)      
2004           continue;
2005         
2006         
2007         AliAODTrack* track = 0;
2008         if(fillNeg)
2009           track = nTrack;
2010         else if(fillPos)
2011           track = pTrack;
2012         else
2013           continue;
2014
2015         Double_t dedx  = track->GetTPCsignal();
2016         Double_t eta  = track->Eta();
2017         Double_t phi  = track->Phi();
2018         Double_t momentum = track->P();
2019
2020         if(track->GetTPCsignalN()<=70)continue;
2021
2022         if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
2023           continue;
2024         
2025         for(Int_t nh = 0; nh < nHists; nh++) {
2026       
2027           if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
2028             continue;
2029
2030           histEV0[nh]->Fill(momentum, dedx);
2031
2032         }
2033         
2034       };
2035         break;
2036
2037
2038       }//end switch
2039     }//end loop over V0s cases
2040     
2041   }//end loop over v0's
2042   
2043   
2044   
2045   
2046 }
2047 Bool_t AliAnalysisTaskQAHighPtDeDx::PhiCut(Double_t pt, Double_t phi, Double_t q, Float_t   mag, TF1* phiCutLow, TF1* phiCutHigh)
2048 {
2049   if(pt < 2.0)
2050     return kTRUE;
2051   
2052   //Double_t phi = track->Phi();
2053   if(mag < 0)    // for negatve polarity field
2054     phi = TMath::TwoPi() - phi;
2055   if(q < 0) // for negatve charge
2056     phi = TMath::TwoPi()-phi;
2057   
2058   phi += TMath::Pi()/18.0; // to center gap in the middle
2059   phi = fmod(phi, TMath::Pi()/9.0);
2060     
2061   if(phi<phiCutHigh->Eval(pt) 
2062      && phi>phiCutLow->Eval(pt))
2063   return kFALSE; // reject track
2064
2065   hPhi->Fill(pt, phi);
2066
2067   return kTRUE;
2068 }