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