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