]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/QATasks/AliAnalysisTaskQAHighPtDeDx.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Completed changes needed because of previous commit
[u/mrichter/AliRoot.git] / PWGLF / QATasks / AliAnalysisTaskQAHighPtDeDx.cxx
CommitLineData
c683985a 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>
62using namespace std;
63
64
65//
66// Responsible:
67// Antonio Ortiz (Lund)
68// Peter Christiansen (Lund)
69//
70
71
72
73
74
75
76
77
78const Double_t AliAnalysisTaskQAHighPtDeDx::fgkClight = 2.99792458e-2;
79Float_t magf = -1;
80TF1* cutLow = new TF1("StandardPhiCutLow", "0.1/x/x+pi/18.0-0.025", 0, 50);
81TF1* cutHigh = new TF1("StandardPhiCutHigh", "0.12/x+pi/18.0+0.035", 0, 50);
82Double_t DeDxMIPMin = 30;
83Double_t DeDxMIPMax = 65;
84const Int_t nHists = 9;
85Float_t centralityGlobal = -10;
86Int_t etaLow[nHists] = {-8, -8, -6, -4, -2, 0, 2, 4, 6};
87Int_t etaHigh[nHists] = { 8, -6, -4, -2, 0, 2, 4, 6, 8};
88
89Int_t nDeltaPiBins = 80;
90Double_t deltaPiLow = 20;
91Double_t deltaPiHigh = 100;
92const Char_t *Pid[7]={"Ch","Pion","Kaon","Proton","Electron","Muon","Oher"};
93ClassImp(AliAnalysisTaskQAHighPtDeDx)
94//_____________________________________________________________________________
95//AliAnalysisTaskQAHighPtDeDx::AliAnalysisTaskQAHighPtDeDx(const char *name):
96AliAnalysisTaskQAHighPtDeDx::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
21cfe1dd 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
c683985a 165}
166
167
168AliAnalysisTaskQAHighPtDeDx::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
240AliAnalysisTaskQAHighPtDeDx::~AliAnalysisTaskQAHighPtDeDx() {
241 //
242 // Destructor
243 //
244
245}
246
247
248
249
250
251
252
253
254//______________________________________________________________________________
255void 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//______________________________________________________________________________
449void 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//________________________________________________________________________
683void 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//________________________________________________________________________
734void 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//_____________________________________________________________________________
776Float_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//_____________________________________________________________________________
789Short_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//_____________________________________________________________________________
819void 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//_____________________________________________________________________________
866void 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//_____________________________________________________________________________
919Short_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//_____________________________________________________________________________
941Short_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//_____________________________________________________________________________
962ULong64_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//____________________________________________________________________
973TParticle* 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//____________________________________________________________________
990Int_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//____________________________________________________________________
1027AliAODMCParticle* 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//____________________________________________________________________
1055TParticle* 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//____________________________________________________________________
1074Int_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//____________________________________________________________________
1114AliAODMCParticle* 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//__________________________________________________________________
1145void 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//__________________________________________________________________
1317void 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
f15c1f69 1326 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*>(AODevent->GetTrack(iT));
1327 if(!aodTrack) AliFatal("Not a standard AOD");
c683985a 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
f15c1f69 1346 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*>(AODevent->GetTrack(iT));
1347 if(!aodTrack) AliFatal("Not a standard AOD");
c683985a 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//__________________________________________________
1504void 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//__________________________________________________________________________
1815void 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}
2047Bool_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}