git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Added PID task, some fixes for coding conventions
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskPWG4PidDetEx.cxx
CommitLineData
de6f8090 1/**************************************************************************
2 * Copyright(c) 1998-1999, 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 <TROOT.h>
17#include <TSystem.h>
18#include <TChain.h>
19#include <TTree.h>
20#include <TFile.h>
21#include <TList.h>
22#include <TStyle.h>
23#include <TCanvas.h>
24#include <TMath.h>
25#include <TH1.h>
26#include <TH2.h>
27#include <TF1.h>
28#include <TProfile.h>
29#include <TLegend.h>
30#include <TDatabasePDG.h>
31#include <TParticle.h>
32
33#include "AliAnalysisManager.h"
34#include "AliAnalysisFilter.h"
35#include "AliESDInputHandler.h"
36#include "AliESDEvent.h"
37#include "AliESDVertex.h"
38#include "AliMCEventHandler.h"
39#include "AliMCEvent.h"
40#include "AliStack.h"
41#include "AliAODInputHandler.h"
42#include "AliAODEvent.h"
43#include "AliAODVertex.h"
44#include "AliAODMCParticle.h"
45#include "AliLog.h"
46
47#include "AliAnalysisTaskPWG4PidDetEx.h"
48
49// STL includes
50#include <iostream>
51using namespace std;
52
53//
54// Analysis class example for PID using detector signals.
55// Works on ESD and AOD; has a flag for MC analysis
56//
57// Alexandru.Dobrin@hep.lu.se
58// Peter.Christiansen@hep.lu.se
59//
60
61ClassImp(AliAnalysisTaskPWG4PidDetEx)
62
63const Double_t AliAnalysisTaskPWG4PidDetEx::fgkCtau = 370; // distance for kaon decay
64const Double_t AliAnalysisTaskPWG4PidDetEx::fgkPionMass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
65const Double_t AliAnalysisTaskPWG4PidDetEx::fgkKaonMass = TDatabasePDG::Instance()->GetParticle(321)->Mass();
66const Double_t AliAnalysisTaskPWG4PidDetEx::fgkProtonMass = TDatabasePDG::Instance()->GetParticle(2212)->Mass();
67
68//_____________________________________________________________________________
69AliAnalysisTaskPWG4PidDetEx::AliAnalysisTaskPWG4PidDetEx():
70 AliAnalysisTaskSE(),
71 fESD(0),
72 fAOD(0),
73 fListOfHists(0),
74 fEtaCut(0.9),
75 fPtCut(0.5),
76 fXbins(20),
77 fXmin(0.),
78 fTOFCutP(2.),
79 fTrackFilter(0x0),
80 fAnalysisMC(kTRUE),
81 fAnalysisType("ESD"),
82 fTriggerMode(kMB2),
83 fEvents(0), fEffTot(0), fEffPID(0), fAccP(0), fAccPt(0), fKaonDecayCorr(0),
84 fdNdPt(0), fMassAll(0),
85 fdNdPtPion(0), fMassPion(0), fdEdxTPCPion(0), fbgTPCPion(0),
86 fdNdPtKaon(0), fMassKaon(0), fdEdxTPCKaon(0), fbgTPCKaon(0),
87 fdNdPtProton(0), fMassProton(0), fdEdxTPCProton(0), fbgTPCProton(0),
88 fdNdPtMC(0), fdNdPtMCPion(0), fdNdPtMCKaon(0), fdNdPtMCProton(0)
89{
90 // Default constructor (should not be used)
91}
92
93//______________________________________________________________________________
94AliAnalysisTaskPWG4PidDetEx::AliAnalysisTaskPWG4PidDetEx(const char *name):
95 AliAnalysisTaskSE(name),
96 fESD(0),
97 fAOD(0),
98 fListOfHists(0),
99 fEtaCut(0.9),
100 fPtCut(0.5),
101 fXbins(20),
102 fXmin(0.),
103 fTOFCutP(2.),
104 fTrackFilter(0x0),
105 fAnalysisMC(kTRUE),
106 fAnalysisType("ESD"),
107 fTriggerMode(kMB2),
108 fEvents(0), fEffTot(0), fEffPID(0), fAccP(0), fAccPt(0), fKaonDecayCorr(0),
109 fdNdPt(0), fMassAll(0),
110 fdNdPtPion(0), fMassPion(0), fdEdxTPCPion(0), fbgTPCPion(0),
111 fdNdPtKaon(0), fMassKaon(0), fdEdxTPCKaon(0), fbgTPCKaon(0),
112 fdNdPtProton(0), fMassProton(0), fdEdxTPCProton(0), fbgTPCProton(0),
113 fdNdPtMC(0), fdNdPtMCPion(0), fdNdPtMCKaon(0), fdNdPtMCProton(0)
114{
115 // Normal constructor
116 // Input slot #0 works with a TChain
117 DefineInput(0, TChain::Class());
118 // Output slot #0 writes into a TList
119 DefineOutput(0, TList::Class());
120
121}
122
123//_____________________________________________________________________________
124AliAnalysisTaskPWG4PidDetEx::~AliAnalysisTaskPWG4PidDetEx()
125{
126 // Destructor
127 // histograms are in the output list and deleted when the output
128 // list is deleted by the TSelector dtor
129 if (fListOfHists) {
130 delete fListOfHists;
131 fListOfHists = 0;
132 }
133}
134
135//______________________________________________________________________________
136void AliAnalysisTaskPWG4PidDetEx::ConnectInputData(Option_t *)
137{
138 // Connect AOD here
139 // Called once
140 if (fDebug > 1) AliInfo("ConnectInputData() \n");
141
142 TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
143 if (!tree) {
144 Printf("ERROR: Could not read chain from input slot 0");
145 }
146 else {
147 if(fAnalysisType == "ESD") {
148 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
149 if (!esdH) {
150 Printf("ERROR: Could not get ESDInputHandler");
151 } else
152 fESD = esdH->GetEvent();
153 }
154 else if(fAnalysisType == "AOD") {
155 AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
156 if (!aodH) {
157 Printf("ERROR: Could not get AODInputHandler");
158 } else
159 fAOD = aodH->GetEvent();
160 }
161 }
162
163}
164
165//______________________________________________________________________________
166void AliAnalysisTaskPWG4PidDetEx::CreateOutputObjects()
167{
168 // Create histograms
169 // Called once
170 if (fDebug > 1) AliInfo("CreateOutPutData() \n");
171
172 OpenFile(0);
173 fListOfHists = new TList();
174
175 fEvents = new TH1I("fEvents","Number of analyzed events; Events; Counts", 1, 0, 1);
176 fListOfHists->Add(fEvents);
177
178 fEffTot = new TH1F ("fEffTot","Efficiency; p_{T} [GeV/c]; Counts", fXbins, fXmin, fTOFCutP);
179 fEffTot->Sumw2();
180 fListOfHists->Add(fEffTot);
181
182 fEffPID = new TH1F ("fEffPID","Efficiency; p_{T} [GeV/c]; Counts", fXbins, fXmin, fTOFCutP);
183 fEffPID->Sumw2();
184 fListOfHists->Add(fEffPID);
185
186 fAccP = new TH1F ("fAccP","Acceptance; p_{T} [GeV/c]; Counts", fXbins, fXmin, fTOFCutP);
187 fAccP->Sumw2();
188 fListOfHists->Add(fAccP);
189
190 fAccPt = new TH1F ("fAccPt","Acceptance; p_{T} [GeV/c]; Counts", fXbins, fXmin, fTOFCutP);
191 fAccPt->Sumw2();
192 fListOfHists->Add(fAccPt);
193
194 fKaonDecayCorr = new TProfile("fKaonDecayCorr","Kaon decay correction vs p_{T}; p_{T} [GeV/c]; Kaon decay correction", fXbins, fXmin, fTOFCutP);
195 fListOfHists->Add(fKaonDecayCorr);
196
197 fdNdPt = new TH1F("fdNdPt","p_{T} distribution; p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", 100 , 0, 10);
198 fdNdPt->Sumw2();
199 fListOfHists->Add(fdNdPt);
200
201 fMassAll = new TH2F("fMassAll","Mass^{2} vs momentum (ToF); p [GeV/c]; M^{2} [GeV^{2}/c^{4}]", 100, -5, 5, 100, -0.2, 1.2);
202 fListOfHists->Add(fMassAll);
203
204 //Pion
205 fdNdPtPion = new TH1F("fdNdPtPion","p_{T} distribution (Pions); p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", fXbins, fXmin, fTOFCutP);
206 fdNdPtPion->Sumw2();
207 fListOfHists->Add(fdNdPtPion);
208
209 fMassPion = new TH2F("fMassPion","Mass squared for pions after cuts vs pmean; p [GeV/c]; Mass^{2} [GeV^{2}/c^{4}]", 100, -5, 5, 100, -0.2, 1.2);
210 fListOfHists->Add(fMassPion);
211
212 fdEdxTPCPion = new TH2F("fdEdxTPCPion","Energy loss vs momentum; p [GeV/c]; dE/dx [a. u.]", 40, -2, 2, 100, 0, 200);
213 fListOfHists->Add(fdEdxTPCPion);
214
215 fbgTPCPion = new TH2F("fbgTPCPion", "Energy loss vs #beta#gamma (Pions);#beta#gamma; dE/dx [a. u.]", 100, 0, 15, 100, 0, 200);
216 fListOfHists->Add(fbgTPCPion);
217
218 //Kaon
219 fdNdPtKaon = new TH1F("fdNdPtKaon","p_{T} distribution (Kaons);p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", fXbins, fXmin, fTOFCutP);
220 fdNdPtKaon->Sumw2();
221 fListOfHists->Add(fdNdPtKaon);
222
223 fMassKaon = new TH2F("fMassKaon","Mass squared for kaons after cuts vs pmean; p [GeV/c]; Mass^{2} [GeV^{2}/c^{4}]", 100, -5, 5, 100, -0.2, 1.2);
224 fListOfHists->Add(fMassKaon);
225
226 fdEdxTPCKaon = new TH2F("fdEdxTPCKaon","Energy loss vs momentum; p [GeV/c]; dE/dx [a. u.]", 40, -2, 2, 100, 0, 200);
227 fListOfHists->Add(fdEdxTPCKaon);
228
229 fbgTPCKaon = new TH2F("fbgTPCKaon", "Energy loss vs #beta#gamma (Kaons);#beta#gamma; dE/dx [a. u.]", 100, 0, 15, 100, 0, 200);
230 fListOfHists->Add(fbgTPCKaon);
231
232 //Proton
233 fdNdPtProton = new TH1F("fdNdPtProton","p_{T} distribution (Protons);p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", fXbins, fXmin, fTOFCutP);
234 fdNdPtProton->Sumw2();
235 fListOfHists->Add(fdNdPtProton);
236
237 fMassProton = new TH2F("fMassProton","Mass squared for protons after cuts vs pmean; p [GeV/c]; Mass^{2} [GeV^{2}/c^{4}]", 100, -5, 5, 100, -0.2, 1.2);
238 fListOfHists->Add(fMassProton);
239
240 fdEdxTPCProton = new TH2F("fdEdxTPCProton","Energy loss vs momentum; p [GeV/c]; dE/dx [a. u.]", 40, -2, 2, 100, 0, 200);
241 fListOfHists->Add(fdEdxTPCProton);
242
243 fbgTPCProton = new TH2F("fbgTPCProton", "Energy loss vs #beta#gamma (Protons);#beta#gamma; dE/dx [a. u.]", 100, 0, 15, 100, 0, 200);
244 fListOfHists->Add(fbgTPCProton);
245
246
247 //MC
248 if(fAnalysisMC){
249 fdNdPtMC = new TH1F("fdNdPtMC","p_{T} distribution;p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", 100, 0, 10);
250 fdNdPtMC->SetLineColor(2);
251 fdNdPtMC->Sumw2();
252 fListOfHists->Add(fdNdPtMC);
253
254 fdNdPtMCPion = new TH1F("fdNdPtMCPion","p_{T} distribution;p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", fXbins, fXmin, fTOFCutP);
255 fdNdPtMCPion->SetLineColor(2);
256 fdNdPtMCPion->Sumw2();
257 fListOfHists->Add(fdNdPtMCPion);
258
259 fdNdPtMCKaon = new TH1F("fdNdPtMCKaon","p_{T} distribution;p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", fXbins, fXmin, fTOFCutP);
260 fdNdPtMCKaon->SetLineColor(2);
261 fdNdPtMCKaon->Sumw2();
262 fListOfHists->Add(fdNdPtMCKaon);
263
264 fdNdPtMCProton = new TH1F("fdNdPtMCProton","p_{T} distribution;p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", fXbins, fXmin, fTOFCutP);
265 fdNdPtMCProton->SetLineColor(2);
266 fdNdPtMCProton->Sumw2();
267 fListOfHists->Add(fdNdPtMCProton);
268 }
269
270}
271
272//______________________________________________________________________________
273void AliAnalysisTaskPWG4PidDetEx::Exec(Option_t *)
274{
275 // Main loop
276 // Called for each event
277 if (fDebug > 1) AliInfo("Exec() \n" );
278
279 if(fAnalysisType == "ESD") {
280 if (!fESD) {
281 Printf("ERROR: fESD not available");
282 return;
283 }
284 if(IsEventTriggered(fESD, fTriggerMode)) {
285 Printf("PWG4 PID ESD analysis");
286 AnalyzeESD(fESD);
287 }//triggered event
288 }//ESD analysis
289
290 else if(fAnalysisType == "AOD") {
291 if (!fAOD) {
292 Printf("ERROR: fAOD not available");
293 return;
294 }
295 if(IsEventTriggered(fAOD, fTriggerMode)) {
296 Printf("PWG4 PID AOD analysis");
297 AnalyzeAOD(fAOD);
298 }//triggered event
299 }//AOD analysis
300
301 // Post output data.
302 PostData(0, fListOfHists);
303}
304
305//________________________________________________________________________
306void AliAnalysisTaskPWG4PidDetEx::AnalyzeESD(AliESDEvent* esd)
307{
308 // Get vertex
309 const AliESDVertex* primaryVertex = esd->GetPrimaryVertex();
310 const Double_t vertexResZ = primaryVertex->GetZRes();
311
312 // Select only events with a reconstructed vertex
313 if(vertexResZ<5.0) {
314 const Int_t nESDTracks = esd->GetNumberOfTracks();
315 for(Int_t i = 0; i < nESDTracks; i++) {
316 AliESDtrack* trackESD = esd->GetTrack(i);
317
318 //Cuts
319 UInt_t selectInfo = 0;
320 if (fTrackFilter) {
321 selectInfo = fTrackFilter->IsSelected(trackESD);
322 if (!selectInfo) continue;
323 }
324
325 if ((trackESD->Pt() < fPtCut) || (TMath::Abs(trackESD->Eta()) > fEtaCut ))
326 continue;
327
328 //Corrections
329 fEffTot->Fill(trackESD->Pt());
330 if (trackESD->GetTOFsignal() !=0)
331 fEffPID->Fill(trackESD->Pt());
332
333 if (trackESD->P() < fTOFCutP) fAccP->Fill(trackESD->Pt());
334 if (trackESD->Pt() < fTOFCutP) fAccPt->Fill(trackESD->Pt());
335
336 //Analysis
337 fdNdPt->Fill(trackESD->Pt(), 1.0/trackESD->Pt());
338
339 //TOF
340 if ((trackESD->GetIntegratedLength() == 0) || (trackESD->GetTOFsignal() == 0))
341 continue;
342
343 fMassAll->Fill(trackESD->Charge()*trackESD->P(), MassSquared(trackESD));
344
345 if ((MassSquared(trackESD) < 0.15) && (MassSquared(trackESD) > -0.15) && (trackESD->P() < fTOFCutP)){
346 fdNdPtPion->Fill(trackESD->Pt(), 1.0/trackESD->Pt());
347 fMassPion->Fill(trackESD->Charge()*trackESD->P(), MassSquared(trackESD));
348 fdEdxTPCPion->Fill(trackESD->Charge()*trackESD->P(),trackESD->GetTPCsignal());
349 fbgTPCPion->Fill(trackESD->P()/fgkPionMass,trackESD->GetTPCsignal());
350 }
351
352 if ((MassSquared(trackESD) > 0.15) && (MassSquared(trackESD) < 0.45) && (trackESD->P() < fTOFCutP)){
353 fdNdPtKaon->Fill(trackESD->Pt(), 1.0/trackESD->Pt());
354 fMassKaon->Fill(trackESD->Charge()*trackESD->P(), MassSquared(trackESD));
355 fdEdxTPCKaon->Fill(trackESD->Charge()*trackESD->P(), trackESD->GetTPCsignal());
356 fbgTPCKaon->Fill(trackESD->P()/fgkKaonMass, trackESD->GetTPCsignal());
357 //Kaon decay correction
358 fKaonDecayCorr->Fill(trackESD->Pt(), TMath::Exp(KaonDecay(trackESD)));
359 }
360
361 if ((MassSquared(trackESD) > 0.75) && (MassSquared(trackESD) < 1.05) && (trackESD->P() < fTOFCutP)){
362 fdNdPtProton->Fill(trackESD->Pt(), 1.0/trackESD->Pt());
363 fMassProton->Fill(trackESD->Charge()*trackESD->P(), MassSquared(trackESD));
364 fdEdxTPCProton->Fill(trackESD->Charge()*trackESD->P(),trackESD->GetTPCsignal());
365 fbgTPCProton->Fill(trackESD->P()/fgkProtonMass,trackESD->GetTPCsignal());
366 }
367 }//ESD track loop
368
369 //MC
370 if(fAnalysisMC){
371 AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
372 if (!mcHandler) {
373 Printf("ERROR: Could not retrieve MC event handler");
374 return;
375 }
376
377 AliMCEvent* mcEvent = mcHandler->MCEvent();
378 if (!mcEvent) {
379 Printf("ERROR: Could not retrieve MC event");
380 return;
381 }
382
383 AliStack* mcStack = mcEvent->Stack();
384 if (!mcStack) {
385 Printf("ERROR: Could not retrieve MC stack");
386 return;
387 }
388
389 const Int_t nTracksMC = mcStack->GetNtrack();
390 for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
391 //Cuts
392 if(!(mcStack->IsPhysicalPrimary(iTracks)))
393 continue;
394
395 TParticle* mcTrack = mcStack->Particle(iTracks);
396
397 Double_t charge = mcTrack->GetPDG()->Charge();
398 if (charge == 0)
399 continue;
400
401 if ((mcTrack->Pt() < fPtCut) || (TMath::Abs(mcTrack->Eta()) > fEtaCut ))
402 continue;
403
404 //Analysis
405 fdNdPtMC->Fill(mcTrack->Pt(), 1.0/mcTrack->Pt());
406
407 if ((mcTrack->GetPdgCode() == 211) || (mcTrack->GetPdgCode() == -211))
408 fdNdPtMCPion->Fill(mcTrack->Pt(),1.0/mcTrack->Pt());
409
410 if ((mcTrack->GetPdgCode() == 321) || (mcTrack->GetPdgCode() == -321))
411 fdNdPtMCKaon->Fill(mcTrack->Pt(),1.0/mcTrack->Pt());
412
413 if ((mcTrack->GetPdgCode() == 2212) || (mcTrack->GetPdgCode() == -2212))
414 fdNdPtMCProton->Fill(mcTrack->Pt(),1.0/mcTrack->Pt());
415 }//MC track loop
416 }//if MC
417 fEvents->Fill(0);
418 }//if vertex
419
420}
421
422//________________________________________________________________________
423void AliAnalysisTaskPWG4PidDetEx::AnalyzeAOD(AliAODEvent* aod)
424{
425 // Get vertex
426 AliAODVertex* primaryVertex = aod->GetPrimaryVertex();
427 const Double_t vertexResZ = TMath::Sqrt(primaryVertex->RotatedCovMatrixZZ());
428
429 // Select only events with a reconstructed vertex
430 if (vertexResZ < 5) {
431 //AOD
432 Int_t nGoodTracks = aod->GetNumberOfTracks();
433
434 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
435 AliAODTrack* trackAOD = aod->GetTrack(iTracks);
436 //Cuts
437 if (!(trackAOD->IsPrimaryCandidate()))
438 continue;
439
440 if ((trackAOD->Pt() < fPtCut) || (TMath::Abs(trackAOD->Eta()) > fEtaCut ))
441 continue;
442
443 //Corrections
444 fEffTot->Fill(trackAOD->Pt());
445 if (trackAOD->GetDetPid()->GetTOFsignal() !=0 )
446 fEffPID->Fill(trackAOD->Pt());
447
448 if (trackAOD->P() < fTOFCutP) fAccP->Fill(trackAOD->Pt());
449 if (trackAOD->Pt() < fTOFCutP) fAccPt->Fill(trackAOD->Pt());
450
451 //Analysis
452 fdNdPt->Fill(trackAOD->Pt(), 1.0/trackAOD->Pt());
453
454 //TOF
455 if ((IntegratedLength(trackAOD) == 0) || (trackAOD->GetDetPid()->GetTOFsignal() == 0))
456 continue;
457
458 fMassAll->Fill(trackAOD->Charge()*trackAOD->P(), MassSquared(trackAOD));
459
460 if ((MassSquared(trackAOD) < 0.15) && (MassSquared(trackAOD) > -0.15) && (trackAOD->P() < fTOFCutP)){
461 fdNdPtPion->Fill(trackAOD->Pt(), 1.0/trackAOD->Pt());
462 fMassPion->Fill(trackAOD->Charge()*trackAOD->P(), MassSquared(trackAOD));
463 fdEdxTPCPion->Fill(trackAOD->Charge()*trackAOD->P(),trackAOD->GetDetPid()->GetTPCsignal());
464 fbgTPCPion->Fill(trackAOD->P()/fgkPionMass,trackAOD->GetDetPid()->GetTPCsignal());
465 }
466
467 if ((MassSquared(trackAOD) > 0.15) && (MassSquared(trackAOD) < 0.45) && (trackAOD->P() < fTOFCutP)){
468 fdNdPtKaon->Fill(trackAOD->Pt(), 1.0/trackAOD->Pt());
469 fMassKaon->Fill(trackAOD->Charge()*trackAOD->P(), MassSquared(trackAOD));
470 fdEdxTPCKaon->Fill(trackAOD->Charge()*trackAOD->P(),trackAOD->GetDetPid()->GetTPCsignal());
471 fbgTPCKaon->Fill(trackAOD->P()/fgkKaonMass,trackAOD->GetDetPid()->GetTPCsignal());
472 //Kaon decay correction
473 fKaonDecayCorr->Fill(trackAOD->Pt(), TMath::Exp(KaonDecay(trackAOD)));
474 }
475
476 if ((MassSquared(trackAOD) > 0.75) && (MassSquared(trackAOD) < 1.05) && (trackAOD->P() < fTOFCutP)){
477 fdNdPtProton->Fill(trackAOD->Pt(), 1.0/trackAOD->Pt());
478 fMassProton->Fill(trackAOD->Charge()*trackAOD->P(), MassSquared(trackAOD));
479 fdEdxTPCProton->Fill(trackAOD->Charge()*trackAOD->P(),trackAOD->GetDetPid()->GetTPCsignal());
480 fbgTPCProton->Fill(trackAOD->P()/fgkProtonMass,trackAOD->GetDetPid()->GetTPCsignal());
481 }
482 }//AOD track loop
483
484 //MC
485 if(fAnalysisMC){
486 TClonesArray* farray = (TClonesArray*)aod->FindListObject("mcparticles");
487 Int_t ntrks = farray->GetEntries();
488 for(Int_t i =0 ; i < ntrks; i++){
489 AliAODMCParticle* trk = (AliAODMCParticle*)farray->At(i);
490 //Cuts
491 if (!(trk->IsPhysicalPrimary()))
492 continue;
493
494 if (trk->Charge() == 0)
495 continue;
496
497 if ((trk->Pt() < fPtCut) || (TMath::Abs(trk->Eta()) > fEtaCut ))
498 continue;
499
500 //Analysis
501 fdNdPtMC->Fill(trk->Pt(), 1.0/trk->Pt());
502
503 if ((trk->GetPdgCode() == 211) || (trk->GetPdgCode() == -211))
504 fdNdPtMCPion->Fill(trk->Pt(),1.0/trk->Pt());
505
506 if ((trk->GetPdgCode() == 321) || (trk->GetPdgCode() == -321))
507 fdNdPtMCKaon->Fill(trk->Pt(),1.0/trk->Pt());
508
509 if ((trk->GetPdgCode() == 2212) || (trk->GetPdgCode() == -2212))
510 fdNdPtMCProton->Fill(trk->Pt(),1.0/trk->Pt());
511 }//MC track loop
512 }//if MC
513 fEvents->Fill(0);
514 }//if vertex resolution
515
516}
517
518//________________________________________________________________________
519Bool_t AliAnalysisTaskPWG4PidDetEx::IsEventTriggered(AliVEvent* ev, TriggerMode trigger)
520{
521 //adapted from PWG2 (AliAnalysisTaskProton)
522
523 ULong64_t triggerMask = ev->GetTriggerMask();
524
525 // definitions from p-p.cfg
526 ULong64_t spdFO = (1 << 14);
527 ULong64_t v0left = (1 << 11);
528 ULong64_t v0right = (1 << 12);
529
530 switch (trigger) {
531 case kMB1:
532 if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
533 return kTRUE;
534 break;
535
536 case kMB2:
537 if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
538 return kTRUE;
539 break;
540
541 case kSPDFASTOR:
542 if (triggerMask & spdFO)
543 return kTRUE;
544 break;
545
546 }//switch
547
548 return kFALSE;
549}
550
551//_____________________________________________________________________________
552Double_t AliAnalysisTaskPWG4PidDetEx::IntegratedLength(AliVTrack* track) const
553{
554 Double_t intTime [5];
555 for (Int_t i = 0; i < 5; i++) intTime[i] = -100.;
556 Double_t timeElectron = 0, intLength = 0;
557
558 AliAODTrack* trackAOD = (AliAODTrack*)track;
559 trackAOD->GetDetPid()->GetIntegratedTimes(intTime);
560 timeElectron = intTime[0];
561 intLength = TMath::C()*timeElectron;
562
563 return intLength;
564}
565
566//_____________________________________________________________________________
567Double_t AliAnalysisTaskPWG4PidDetEx::MassSquared(AliVTrack* track) const
568{
569 Double_t beta = -10, mass = -10;
570
571 if(fAnalysisType == "ESD"){
572 AliESDtrack* trackESD = (AliESDtrack*)track;
573 beta = trackESD->GetIntegratedLength()/trackESD->GetTOFsignal()/TMath::C();
574 mass = trackESD->P()*trackESD->P()*(1./(beta*beta) - 1.0);
575 }
576
577 if(fAnalysisType == "AOD"){
578 AliAODTrack* trackAOD = (AliAODTrack*)track;
579 beta =IntegratedLength(trackAOD)/trackAOD->GetDetPid()->GetTOFsignal()/TMath::C();
580 mass = trackAOD->P()*trackAOD->P()*(1./(beta*beta) - 1.0);
581 }
582
583 return mass;
584}
585
586//_____________________________________________________________________________
587Double_t AliAnalysisTaskPWG4PidDetEx::KaonDecay(AliVTrack* track) const
588{
589 Double_t decay = -10;
590
591 if(fAnalysisType == "ESD"){
592 AliESDtrack* trackESD = (AliESDtrack*)track;
593 decay = trackESD->GetIntegratedLength()*fgkKaonMass/fgkCtau/trackESD->P();
594 }
595
596 if (fAnalysisType == "AOD"){
597 AliAODTrack* trackAOD = (AliAODTrack*)track;
598 decay = IntegratedLength(trackAOD)*fgkKaonMass/fgkCtau/trackAOD->P();
599 }
600
601 return decay;
602}
603
604//_____________________________________________________________________________
605void AliAnalysisTaskPWG4PidDetEx::Terminate(Option_t *)
606{
607 // Terminate loop
608 if (fDebug > 1) Printf("Terminate()");
609
610 //
611 // The followig code has now been moved to drawPid.C
612 //
613
614// fListOfHists = dynamic_cast<TList*> (GetOutputData(0));
615// if (!fListOfHists) {
616// Printf("ERROR: fListOfHists not available");
617// return;
618// }
619
620// TH1I* hevents = dynamic_cast<TH1I*> (fListOfHists->FindObject("fEvents"));
621// Int_t nEvents = (Int_t) hevents->GetBinContent(1);
622
623// const Float_t normalization = 2.0*fEtaCut*2.0*TMath::Pi();
624
625// //-----------------
626// TCanvas* c1 = new TCanvas("c1", "c1");
627// c1->cd();
628
629// TH2F* hMassAll = dynamic_cast<TH2F*> (fListOfHists->FindObject("fMassAll"));
630// hMassAll->SetLineColor(1);
631// hMassAll->SetMarkerColor(1);
632// hMassAll->DrawCopy();
633
634// TH2F* hMassPion = dynamic_cast<TH2F*> (fListOfHists->FindObject("fMassPion"));
635// hMassPion->SetLineColor(2);
636// hMassPion->SetMarkerColor(2);
637// hMassPion->SetMarkerStyle(7);
638// hMassPion->DrawCopy("same");
639
640// TH2F* hMassKaon = dynamic_cast<TH2F*> (fListOfHists->FindObject("fMassKaon"));
641// hMassKaon->SetLineColor(3);
642// hMassKaon->SetMarkerColor(3);
643// hMassKaon->SetMarkerStyle(7);
644// hMassKaon->DrawCopy("same");
645
646// TH2F* hMassProton = dynamic_cast<TH2F*> (fListOfHists->FindObject("fMassProton"));
647// hMassProton->SetLineColor(4);
648// hMassProton->SetMarkerColor(4);
649// hMassProton->SetMarkerStyle(7);
650// hMassProton->DrawCopy("same");
651
652// TLegend* legend1 = new TLegend(0.8, 0.8, 1, 1);
653// legend1->SetBorderSize(0);
654// legend1->SetFillColor(0);
655// legend1->AddEntry(hMassAll, "All","LP");
656// legend1->AddEntry(hMassPion, "Pions","LP");
657// legend1->AddEntry(hMassKaon, "Kaons","LP");
658// legend1->AddEntry(hMassProton, "Protons","LP");
659// legend1->Draw();
660
661// c1->Update();
662
663// //-----------------
664// TCanvas* c2 = new TCanvas("c2", "c2");
665// c2->cd();
666
667// TH2F* hdEdxTPCPion = dynamic_cast<TH2F*> (fListOfHists->FindObject("fdEdxTPCPion"));
668// hdEdxTPCPion->SetTitle("dE/dx vs p (TPC)");
669// hdEdxTPCPion->SetLineColor(2);
670// hdEdxTPCPion->SetMarkerColor(2);
671// hdEdxTPCPion->SetMarkerStyle(7);
672// hdEdxTPCPion->DrawCopy();
673
674// TH2F* hdEdxTPCKaon = dynamic_cast<TH2F*> (fListOfHists->FindObject("fdEdxTPCKaon"));
675// hdEdxTPCKaon->SetLineColor(3);
676// hdEdxTPCKaon->SetMarkerColor(3);
677// hdEdxTPCKaon->SetMarkerStyle(7);
678// hdEdxTPCKaon->DrawCopy("same");
679
680// TH2F* hdEdxTPCProton = dynamic_cast<TH2F*> (fListOfHists->FindObject("fdEdxTPCProton"));
681// hdEdxTPCProton->SetLineColor(4);
682// hdEdxTPCProton->SetMarkerColor(4);
683// hdEdxTPCProton->SetMarkerStyle(7);
684// hdEdxTPCProton->DrawCopy("same");
685
686// TLegend* legend2 = new TLegend(0.66, 0.66, 0.88, 0.88);
687// legend2->SetBorderSize(0);
688// legend2->SetFillColor(0);
689// legend2->AddEntry(hdEdxTPCPion, "Pions","LP");
690// legend2->AddEntry(hdEdxTPCKaon, "Kaons","LP");
691// legend2->AddEntry(hdEdxTPCProton, "Protons","LP");
692// legend2->Draw();
693
694// c2->Update();
695
696// //-----------------
697// TCanvas* c3 = new TCanvas("c3", "c3");
698// c3->cd();
699
700// TH2F* hdEdxTPCbgPion = dynamic_cast<TH2F*> (fListOfHists->FindObject("fbgTPCPion"));
701// hdEdxTPCbgPion->SetTitle("dE/dx vs #beta#gamma (TPC)");
702// hdEdxTPCbgPion->SetLineColor(2);
703// hdEdxTPCbgPion->SetMarkerColor(2);
704// hdEdxTPCbgPion->SetMarkerStyle(7);
705// hdEdxTPCbgPion->DrawCopy();
706
707// TH2F* hdEdxTPCbgKaon = dynamic_cast<TH2F*> (fListOfHists->FindObject("fbgTPCKaon"));
708// hdEdxTPCbgKaon->SetLineColor(3);
709// hdEdxTPCbgKaon->SetMarkerColor(3);
710// hdEdxTPCbgKaon->SetMarkerStyle(7);
711// hdEdxTPCbgKaon->DrawCopy("same");
712
713// TH2F* hdEdxTPCbgProton = dynamic_cast<TH2F*> (fListOfHists->FindObject("fbgTPCProton"));
714// hdEdxTPCbgProton->SetLineColor(4);
715// hdEdxTPCbgProton->SetMarkerColor(4);
716// hdEdxTPCbgProton->SetMarkerStyle(7);
717// hdEdxTPCbgProton->DrawCopy("same");
718
719// TLegend* legend3 = new TLegend(0.66, 0.66, 0.88, 0.88);
720// legend3->SetBorderSize(0);
721// legend3->SetFillColor(0);
722// legend3->AddEntry(hdEdxTPCbgPion, "Pions","LP");
723// legend3->AddEntry(hdEdxTPCbgKaon, "Kaons","LP");
724// legend3->AddEntry(hdEdxTPCbgProton, "Protons","LP");
725// legend3->Draw();
726
727// c3->Update();
728
729// //-----------------
730// TCanvas* c4 = new TCanvas("c4", "c4", 100, 100, 500, 900);
731// c4->Divide(1,3);
732
733// c4->cd(1);
734// TH1F* hAccepatncePt = dynamic_cast<TH1F*> (fListOfHists->FindObject("fAccPt"));
735// TH1F* hAcceptanceP = dynamic_cast<TH1F*> (fListOfHists->FindObject("fAccP"));
736// hAccepatncePt->Divide(hAcceptanceP);
737// hAccepatncePt->SetTitle("Acceptance correction");
738// hAccepatncePt->GetYaxis()->SetTitle("Acceptance correction");
739// hAccepatncePt->DrawCopy();
740
741// c4->cd(2);
742// TH1F* hEfficiencyPID = dynamic_cast<TH1F*> (fListOfHists->FindObject("fEffPID"));
743// TH1F* hEfficiencyTot = dynamic_cast<TH1F*> (fListOfHists->FindObject("fEffTot"));
744// hEfficiencyPID->Divide(hEfficiencyTot);
745// hEfficiencyPID->SetTitle("Efficiency correction");
746// hEfficiencyPID->GetYaxis()->SetTitle("Efficiency correction");
747// hEfficiencyPID->DrawCopy();
748
749// c4->cd(3);
750// TProfile* hKDecayCorr = dynamic_cast<TProfile*> (fListOfHists->FindObject("fKaonDecayCorr"));
751// hKDecayCorr->SetTitle("Kaon decay correction");
752// hKDecayCorr->GetYaxis()->SetTitle("Kaon decay correction");
753// hKDecayCorr->GetXaxis()->SetTitle(" p_{T} [GeV/c]");
754// hKDecayCorr->DrawCopy();
755
756// c4->Update();
757
758// //---------------------
759// TCanvas* c5 = new TCanvas("c5", "c5", 100, 100, 600, 900);
760// c5->Divide(1,2);
761
762// c5->cd(1);
763// TH1F* hPtPionNoCorr = dynamic_cast<TH1F*> (fListOfHists->FindObject("fdNdPtPion"));
764// hPtPionNoCorr->Scale(1.0/nEvents/normalization/hPtPionNoCorr->GetBinWidth(1));
765// hPtPionNoCorr->SetTitle("p_{T} distribution (no corrections)");
766// hPtPionNoCorr->SetLineColor(2);
767// hPtPionNoCorr->GetYaxis()->SetRangeUser(1E-3,1E1);
768// hPtPionNoCorr->DrawCopy();
769
770// TH1F* hPtKaonNoCorr = dynamic_cast<TH1F*> (fListOfHists->FindObject("fdNdPtKaon"));
771// hPtKaonNoCorr->Scale(1.0/nEvents/normalization/hPtKaonNoCorr->GetBinWidth(1));
772// hPtKaonNoCorr->SetLineColor(3);
773// hPtKaonNoCorr->GetYaxis()->SetRangeUser(1E-3,1E1);
774// hPtKaonNoCorr->DrawCopy("same");
775
776// TH1F* hPtProtonNoCorr = dynamic_cast<TH1F*> (fListOfHists->FindObject("fdNdPtProton"));
777// hPtProtonNoCorr->Scale(1.0/nEvents/normalization/hPtProtonNoCorr->GetBinWidth(1));
778// hPtProtonNoCorr->SetLineColor(4);
779// hPtProtonNoCorr->GetYaxis()->SetRangeUser(1E-3,1E1);
780// hPtProtonNoCorr->DrawCopy("same");
781
782// TH1F* hPt = dynamic_cast<TH1F*> (fListOfHists->FindObject("fdNdPt"));
783// hPt->Scale(1.0/nEvents/normalization/hPt->GetBinWidth(1));
784// hPt->GetYaxis()->SetRangeUser(1E-3,1E1);
785// hPt->DrawCopy("same");
786
787// TLegend* legend4 = new TLegend(0.63, 0.63, 0.88, 0.88);
788// legend4->SetBorderSize(0);
789// legend4->SetFillColor(0);
790// legend4->AddEntry(hPt, "p_{T} dist (all)","L");
791// legend4->AddEntry(hPtPionNoCorr, "p_{T} dist (Pions)","L");
792// legend4->AddEntry(hPtKaonNoCorr, "p_{T} dist (Kaons)","L");
793// legend4->AddEntry(hPtProtonNoCorr, "p_{T} dist (Protons)","L");
794// legend4->Draw();
795
796// gPad->SetLogy();
797
798
799// c5->cd(2);
800// TH1F* hPtPionCorr = static_cast<TH1F*>(hPtPionNoCorr->Clone());
801// hPtPionCorr->SetTitle("p_{T} distribution (with corrections)");
802// hPtPionCorr->Multiply(hAccepatncePt);
803// hPtPionCorr->Divide(hEfficiencyPID);
804// hPtPionCorr->GetYaxis()->SetRangeUser(1E-2,1E1);
805// hPtPionCorr->DrawCopy();
806
807// TH1F* hPtKaonCorr = static_cast<TH1F*>(hPtKaonNoCorr->Clone());
808// hPtKaonCorr->Multiply(hAccepatncePt);
809// hPtKaonCorr->Divide(hEfficiencyPID);
810// hPtKaonCorr->Multiply(hKDecayCorr);
811// hPtKaonCorr->GetYaxis()->SetRangeUser(1E-2,1E1);
812// hPtKaonCorr->DrawCopy("same");
813
814// TH1F* hPtProtonCorr = static_cast<TH1F*>(hPtProtonNoCorr->Clone());
815// hPtProtonCorr->Multiply(hAccepatncePt);
816// hPtProtonCorr->Divide(hEfficiencyPID);
817// hPtProtonCorr->GetYaxis()->SetRangeUser(1E-2,1E1);
818// hPtProtonCorr->DrawCopy("same");
819
820// hPt->GetYaxis()->SetRangeUser(1E-2,1E1);
821// hPt->DrawCopy("same");
822
823// TLegend* legend5 = new TLegend(0.63, 0.63, 0.88, 0.88);
824// legend5->SetBorderSize(0);
825// legend5->SetFillColor(0);
826// legend5->AddEntry(hPt, "p_{T} dist (all)","L");
827// legend5->AddEntry(hPtPionCorr, "p_{T} dist (Pions)","L");
828// legend5->AddEntry(hPtKaonCorr, "p_{T} dist (Kaons)","L");
829// legend5->AddEntry(hPtProtonCorr, "p_{T} dist (Protons)","L");
830// legend5->Draw();
831
832// gPad->SetLogy();
833
834// c5->Update();
835
836// //-----------------
837// if (fAnalysisMC){
838// TCanvas* c6 = new TCanvas("c6", "c6", 100, 100, 1200, 800);
839// c6->Divide(2,2);
840
841// c6->cd(1);
842// TH1F* hPt_clone = static_cast<TH1F*>(hPt->Clone());
843// hPt_clone->GetYaxis()->SetRangeUser(1E-5,1E1);
844// hPt_clone->SetTitle("p_{T} distribution (all)");
845// hPt_clone->SetLineColor(1);
846// hPt_clone->DrawCopy();
847
848// TH1F* hPtMC = dynamic_cast<TH1F*> (fListOfHists->FindObject("fdNdPtMC"));
849// hPtMC->Scale(1.0/nEvents/normalization/hPtMC->GetBinWidth(1));
850// hPtMC->GetYaxis()->SetRangeUser(1E-5,1E1);
851// hPtMC->DrawCopy("same");
852
853// TLegend* legend6 = new TLegend(0.57, 0.57, 0.87, 0.87);
854// legend6->SetBorderSize(0);
855// legend6->SetFillColor(0);
856// legend6->AddEntry(hPt_clone, "p_{T} dist (Rec)", "L");
857// legend6->AddEntry(hPtMC, "p_{T} dist (MC)", "L");
858// legend6->Draw();
859
860// gPad->SetLogy();
861
862// c6->cd(2);
863// TH1F* hPtPion_clone = static_cast<TH1F*>(hPtPionCorr->Clone());
864// hPtPion_clone->GetYaxis()->SetRangeUser(1E-2,1E1);
865// hPtPion_clone->SetTitle("p_{T} distribution (Pions)");
866// hPtPion_clone->SetLineColor(1);
867// hPtPion_clone->DrawCopy();
868
869// TH1F* hPtMCPion = dynamic_cast<TH1F*> (fListOfHists->FindObject("fdNdPtMCPion"));
870// hPtMCPion->Scale(1.0/nEvents/normalization/hPtMCPion->GetBinWidth(1));
871// hPtMCPion->GetYaxis()->SetRangeUser(1E-2,1E1);
872// hPtMCPion->DrawCopy("same");
873
874// TLegend* legend7 = new TLegend(0.57, 0.57, 0.87, 0.87);
875// legend7->SetBorderSize(0);
876// legend7->SetFillColor(0);
877// legend7->AddEntry(hPtPion_clone, "p_{T} dist (Rec)", "L");
878// legend7->AddEntry(hPtMCPion, "p_{T} dist (MC)", "L");
879// legend7->Draw();
880
881// gPad->SetLogy();
882
883// c6->cd(3);
884// TH1F* hPtKaon_clone = static_cast<TH1F*>(hPtKaonCorr->Clone());
885// hPtKaon_clone->GetYaxis()->SetRangeUser(1E-2,1E0);
886// hPtKaon_clone->SetLineColor(1);
887// hPtKaon_clone->DrawCopy();
888
889// TH1F* hPtMCKaon = dynamic_cast<TH1F*> (fListOfHists->FindObject("fdNdPtMCKaon"));
890// hPtMCKaon->Scale(1.0/nEvents/normalization/hPtMCKaon->GetBinWidth(1));
891// hPtMCKaon->GetYaxis()->SetRangeUser(1E-2,1E0);
892// hPtMCKaon->DrawCopy("same");
893
894// TLegend* legend8 = new TLegend(0.57, 0.57, 0.87, 0.87);
895// legend8->SetBorderSize(0);
896// legend8->SetFillColor(0);
897// legend8->AddEntry(hPtKaon_clone, "p_{T} dist (Rec)", "L");
898// legend8->AddEntry(hPtMCKaon, "p_{T} dist (MC)", "L");
899// legend8->Draw();
900
901// gPad->SetLogy();
902
903// c6->cd(4);
904// TH1F* hPtProton_clone = static_cast<TH1F*>(hPtProtonCorr->Clone());
905// hPtProton_clone->GetYaxis()->SetRangeUser(1E-2,1E-1);
906// hPtProton_clone->SetLineColor(1);
907// hPtProton_clone->DrawCopy();
908
909// TH1F* hPtMCProton = dynamic_cast<TH1F*> (fListOfHists->FindObject("fdNdPtMCProton"));
910// hPtMCProton->Scale(1.0/nEvents/normalization/hPtMCProton->GetBinWidth(1));
911// hPtMCProton->GetYaxis()->SetRangeUser(1E-2,1E-1);
912// hPtMCProton->DrawCopy("same");
913
914// TLegend* legend9 = new TLegend(0.2, 0.25, 0.5, 0.55);
915// legend9->SetBorderSize(0);
916// legend9->SetFillColor(0);
917// legend9->AddEntry(hPtProton_clone, "p_{T} dist (Rec)", "L");
918// legend9->AddEntry(hPtMCProton, "p_{T} dist (MC)", "L");
919// legend9->Draw();
920
921// gPad->SetLogy();
922
923// c6->Update();
924// }
925
926}