]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/Correlations/DPhi/DiHadronPID/AliAnalysisTaskDiHadronPID.cxx
modification on the usage of some pid-flags (A.Rossi)
[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / DiHadronPID / AliAnalysisTaskDiHadronPID.cxx
CommitLineData
97724bd1 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// -----------------------------------------------------------------------
17// This analysis task fills histograms with PID information of tracks
18// associated to a high p_T trigger.
19// -----------------------------------------------------------------------
20// Author: Misha Veldhoen (misha.veldhoen@cern.ch)
6788af99 21
22#include <iostream>
23
24// Basic Includes
25#include "TH1F.h"
26#include "TH2F.h"
27#include "TH3F.h"
28#include "THn.h"
29#include "TFile.h"
30#include "TChain.h"
31#include "TObject.h"
32#include "TObjArray.h"
33
34// Manager/ Handler
35#include "AliAnalysisManager.h"
36#include "AliInputEventHandler.h"
37
38// Event pool includes.
39#include "AliEventPoolManager.h"
40
41// PID includes.
42#include "AliPIDResponse.h"
43
44// AOD includes.
45#include "AliAODEvent.h"
46#include "AliAODTrack.h"
47#include "AliAODHandler.h"
48#include "AliAODVertex.h"
49#include "AliAODInputHandler.h"
50#include "AliAODMCParticle.h"
51#include "AliAODMCHeader.h"
52
53// Additional includes.
54#include "AliTrackDiHadronPID.h"
55#include "AliAODTrackCutsDiHadronPID.h"
56#include "AliAODEventCutsDiHadronPID.h"
57#include "AliHistToolsDiHadronPID.h"
5c01a71f 58#include "AliFunctionsDiHadronPID.h"
6788af99 59
60// AnalysisTask headers.
61#include "AliAnalysisTaskSE.h"
62#include "AliAnalysisTaskDiHadronPID.h"
63
64using namespace std;
65
66ClassImp(AliAnalysisTaskDiHadronPID);
67
97724bd1 68// -----------------------------------------------------------------------
6788af99 69AliAnalysisTaskDiHadronPID::AliAnalysisTaskDiHadronPID():
70 AliAnalysisTaskSE(),
71 fPIDResponse(0x0),
72 fEventCuts(0x0),
73 fTrackCutsTrigger(0x0),
74 fTrackCutsAssociated(0x0),
75 fPoolMgr(0x0),
76 fTriggerTracks(0x0),
77 fAssociatedTracks(0x0),
78 fCurrentAODEvent(0x0),
79 fOutputList(0x0),
07d62e30 80 fPtSpectrumTOFbins(0x0),
81 fCorrelationsTOFbins(0x0),
82 fMixedEventsTOFbins(0x0),
83 fPtSpectrumTOFTPCbins(0x0),
84 fCorrelationsTOFTPCbins(0x0),
85 fMixedEventsTOFTPCbins(0x0),
6788af99 86 fTOFhistos(0x0),
07d62e30 87 fTOFmismatch(0x0),
5c01a71f 88 fTOFPtAxis(0x0),
89 fTOFTPChistos(0x0),
07d62e30 90 fTOFTPCmismatch(0x0),
5c01a71f 91 fTOFTPCPtAxis(0x0),
6788af99 92 fNDEtaBins(32),
93 fNDPhiBins(32),
94 fMinNEventsForMixing(5),
95 fPoolTrackDepth(2000),
96 fPoolSize(1000),
50dfda71 97 fMixEvents(kTRUE),
98 fMixTriggers(kFALSE),
07d62e30 99 fCalculateMismatch(kTRUE),
6788af99 100 fT0Fill(0x0),
101 fLvsEta(0x0),
5c01a71f 102 fLvsEtaProjections(0x0),
103 fMakeTOFcorrelations(kTRUE),
0770be27 104 fMakeTOFTPCcorrelationsPi(kFALSE),
105 fMakeTOFTPCcorrelationsKa(kFALSE),
106 fMakeTOFTPCcorrelationsPr(kFALSE),
107 fTOFIntervalFactorTOFTPC(1.)
6788af99 108
109{
110
111 //
112 // Default Constructor.
113 //
114
115 if (fDebug > 0) {AliInfo("AliAnalysisTaskDiHadronPID Default Constructor.");}
116
6788af99 117}
118
97724bd1 119// -----------------------------------------------------------------------
6788af99 120AliAnalysisTaskDiHadronPID::AliAnalysisTaskDiHadronPID(const char* name):
121 AliAnalysisTaskSE(name),
122 fPIDResponse(0x0),
123 fEventCuts(0x0),
124 fTrackCutsTrigger(0x0),
125 fTrackCutsAssociated(0x0),
126 fPoolMgr(0x0),
127 fTriggerTracks(0x0),
128 fAssociatedTracks(0x0),
129 fCurrentAODEvent(0x0),
130 fOutputList(0x0),
07d62e30 131 fPtSpectrumTOFbins(0x0),
132 fCorrelationsTOFbins(0x0),
133 fMixedEventsTOFbins(0x0),
134 fPtSpectrumTOFTPCbins(0x0),
135 fCorrelationsTOFTPCbins(0x0),
136 fMixedEventsTOFTPCbins(0x0),
5c01a71f 137 fTOFhistos(0x0),
07d62e30 138 fTOFmismatch(0x0),
5c01a71f 139 fTOFPtAxis(0x0),
140 fTOFTPChistos(0x0),
07d62e30 141 fTOFTPCmismatch(0x0),
142 fTOFTPCPtAxis(0x0),
6788af99 143 fNDEtaBins(32),
144 fNDPhiBins(32),
145 fMinNEventsForMixing(5),
146 fPoolTrackDepth(2000),
147 fPoolSize(1000),
50dfda71 148 fMixEvents(kTRUE),
149 fMixTriggers(kFALSE),
07d62e30 150 fCalculateMismatch(kTRUE),
6788af99 151 fT0Fill(0x0),
152 fLvsEta(0x0),
5c01a71f 153 fLvsEtaProjections(0x0),
154 fMakeTOFcorrelations(kTRUE),
0770be27 155 fMakeTOFTPCcorrelationsPi(kFALSE),
156 fMakeTOFTPCcorrelationsKa(kFALSE),
157 fMakeTOFTPCcorrelationsPr(kFALSE),
158 fTOFIntervalFactorTOFTPC(1.)
6788af99 159
160{
161
162 //
163 // Named Constructor.
164 //
165
50dfda71 166 if (fDebug > 0) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
6788af99 167
6788af99 168 DefineInput(0,TChain::Class());
169 DefineOutput(1,TList::Class());
170
171}
172
97724bd1 173// -----------------------------------------------------------------------
6788af99 174AliAnalysisTaskDiHadronPID::~AliAnalysisTaskDiHadronPID() {;
175
176 //
177 // Destructor.
178 //
179
50dfda71 180 if (fDebug > 0) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
6788af99 181
5c01a71f 182 if (fPoolMgr) {delete fPoolMgr; fPoolMgr = 0x0;}
183 if (fOutputList) {delete fOutputList; fOutputList = 0x0;}
184
6788af99 185}
186
97724bd1 187// -----------------------------------------------------------------------
6788af99 188void AliAnalysisTaskDiHadronPID::UserCreateOutputObjects() {
189
190 //
191 // Create Output objects.
192 //
193
50dfda71 194 if (fDebug > 0) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
6788af99 195
196 // --- BEGIN: Initialization on the worker nodes ---
197 AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
198 AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*> (manager->GetInputEventHandler());
199
200 // Getting the pointer to the PID response object.
201 fPIDResponse = inputHandler->GetPIDResponse();
202
07d62e30 203 // For now we don't bin in multiplicity for pp.
ae8330d1 204 TArrayD* centralityBins = 0x0;
07d62e30 205 if (fEventCuts->GetIsPbPb()) {
07d62e30 206 Double_t tmp[] = {0., 1., 2., 3., 4., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.1 };
ae8330d1 207 centralityBins = new TArrayD(15, tmp);
07d62e30 208 } else {
07d62e30 209 Double_t tmp[] = {0.,1.};
ae8330d1 210 centralityBins = new TArrayD(2, tmp);
07d62e30 211 }
6788af99 212
213 Int_t nZvtxBins = 7;
5c01a71f 214 Double_t vertexBins[] = {-7., -5., -3., -1., 1., 3., 5., 7.};
6788af99 215
ae8330d1 216 fPoolMgr = new AliEventPoolManager(fPoolSize, fPoolTrackDepth, centralityBins->GetSize(), centralityBins->GetArray(), nZvtxBins, (Double_t*) vertexBins);
217
218 delete centralityBins;
6788af99 219 // --- END ---
220
221 // Create the output list.
222 fOutputList = new TList();
223 fOutputList->SetOwner(kTRUE);
224
225 // Creating all requested histograms locally.
226 fEventCuts->CreateHistos();
227 fOutputList->Add(fEventCuts);
228
229 fTrackCutsTrigger->CreateHistos();
230 fOutputList->Add(fTrackCutsTrigger);
231
232 fTrackCutsAssociated->CreateHistos();
233 fOutputList->Add(fTrackCutsAssociated);
234
5c01a71f 235 TString speciesname[] = {"Pion","Kaon","Proton"};
6788af99 236
5c01a71f 237 // Create TOF correlations histograms (DPhi,DEta,TOF).
238 if (fMakeTOFcorrelations) {
6788af99 239
07d62e30 240 // Get the pT axis for the TOF PID correlations.
241 Double_t* ptaxis = fTrackCutsAssociated->GetPtAxisPID();
242 Int_t nptbins = fTrackCutsAssociated->GetNPtBinsPID();
243
244 // Create Pt spectrum histogram.
245 fPtSpectrumTOFbins = new TH1F("fPtSpectrumTOFbins","p_{T} Spectrum;p_{T} (GeV/c);Count",nptbins,ptaxis);
246 fOutputList->Add(fPtSpectrumTOFbins);
247
248 // Create unidentified correlations histogram.
249 fCorrelationsTOFbins = AliHistToolsDiHadronPID::MakeHist3D("fCorrelationsTOFbins","Correlations;#Delta#phi;#Delta#eta;p_{T} (GeV/c)",
250 fNDPhiBins,-TMath::Pi()/2.,3.*TMath::Pi()/2.,
251 fNDEtaBins,-1.6,1.6,
252 nptbins, ptaxis);
253 fOutputList->Add(fCorrelationsTOFbins);
254
255 // Create unidentified mixed events histogram.
256 fMixedEventsTOFbins = AliHistToolsDiHadronPID::MakeHist3D("fMixedEventsTOFbins","Mixed Events;#Delta#phi;#Delta#eta;p_{T} (GeV/c)",
257 fNDPhiBins,-TMath::Pi()/2.,3.*TMath::Pi()/2.,
258 fNDEtaBins,-1.6,1.6,
259 nptbins, ptaxis);
260 fOutputList->Add(fMixedEventsTOFbins);
261
262 // Create TOFPtaxis.
263 fTOFPtAxis = new TAxis(nptbins, ptaxis);
264 fTOFPtAxis->SetName("fTOFPtAxis");
265 fTOFPtAxis->SetTitle("p_{T} GeV/c");
266
267 // Create PID histograms.
5c01a71f 268 fTOFhistos = new TObjArray(3);
269 fTOFhistos->SetOwner(kTRUE);
270 fTOFhistos->SetName("CorrelationsTOF");
6788af99 271
07d62e30 272 if (fCalculateMismatch) {
273 fTOFmismatch = new TObjArray(3);
274 fTOFmismatch->SetOwner(kTRUE);
275 fTOFmismatch->SetName("MismatchTOF");
276 }
277
6788af99 278 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
279
07d62e30 280 TObjArray* TOFhistosTmp = new TObjArray(fTOFPtAxis->GetNbins());
281 TOFhistosTmp->SetOwner(kTRUE);
282 TOFhistosTmp->SetName(speciesname[iSpecies].Data());
283
284 TObjArray* TOFmismatchTmp = 0x0;
285 if (fCalculateMismatch) {
286 TOFmismatchTmp = new TObjArray(fTOFPtAxis->GetNbins());
287 TOFmismatchTmp->SetOwner(kTRUE);
288 TOFmismatchTmp->SetName(speciesname[iSpecies].Data());
289 }
5c01a71f 290
291 for (Int_t iBinPt = 1; iBinPt < (fTOFPtAxis->GetNbins() + 1); iBinPt++) {
292
293 Int_t iPtClass = fTrackCutsAssociated->GetPtClass(iBinPt);
ae8330d1 294 if (iPtClass == -1) {AliFatal("Not valid pT class.");}
5c01a71f 295
296 Int_t NBinsTOF = fTrackCutsAssociated->GetNTOFbins(iPtClass,iSpecies);
297 Double_t TOFmin = fTrackCutsAssociated->GetTOFmin(iPtClass,iSpecies);
298 Double_t TOFmax = fTrackCutsAssociated->GetTOFmax(iPtClass,iSpecies);
299
07d62e30 300 //cout << "ptbin: "<< iBinPt << " class: " << iPtClass << " TOFBins: " << NBinsTOF << " min: " << TOFmin << " max: " << TOFmax << endl;
6788af99 301
07d62e30 302 // Correlation histogram.
5c01a71f 303 TH3F* htmp = new TH3F(Form("fCorrelationsTOF_%i",iBinPt),
304 Form("%5.3f < p_{T} < %5.3f; #Delta#phi; #Delta#eta; t_{TOF} (ps)", fTOFPtAxis->GetBinLowEdge(iBinPt), fTOFPtAxis->GetBinUpEdge(iBinPt)),
305 fNDPhiBins, -TMath::Pi()/2., 3.*TMath::Pi()/2.,
306 fNDEtaBins, -1.6, 1.6, NBinsTOF, TOFmin, TOFmax);
307 htmp->SetDirectory(0);
6788af99 308
07d62e30 309 TOFhistosTmp->Add(htmp);
5c01a71f 310
07d62e30 311 if (fCalculateMismatch) {
312 // Mismatch histogram.
313 TH1F* htmp2 = new TH1F(Form("fMismatchTOF_%i",iBinPt),
314 Form("%5.3f < p_{T} < %5.3f; t_{TOF} (ps)", fTOFPtAxis->GetBinLowEdge(iBinPt), fTOFPtAxis->GetBinUpEdge(iBinPt)),
315 NBinsTOF, TOFmin, TOFmax);
316 htmp2->SetDirectory(0);
317
318 TOFmismatchTmp->Add(htmp2);
319 }
5c01a71f 320 }
321
07d62e30 322 fTOFhistos->Add(TOFhistosTmp);
323 if (fCalculateMismatch) {fTOFmismatch->Add(TOFmismatchTmp);}
6788af99 324
325 }
5c01a71f 326
327 fOutputList->Add(fTOFhistos);
07d62e30 328 if (fCalculateMismatch) {fOutputList->Add(fTOFmismatch);}
5c01a71f 329
6788af99 330 }
331
5c01a71f 332 // Create TOF/TPC correlation histograms. (DPhi,DEta,TOF,TPC).
0770be27 333 if (fMakeTOFTPCcorrelationsPi || fMakeTOFTPCcorrelationsKa || fMakeTOFTPCcorrelationsPr) {
5c01a71f 334
335 Double_t ptarrayTOFTPC[16] = {2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6,
336 2.8, 3.0, 3.2, 3.4, 3.6, 3.8,
337 4.2, 4.6, 5.0};
07d62e30 338 const Int_t nptbins = 15;
339 fTOFTPCPtAxis = new TAxis(nptbins, ptarrayTOFTPC);
5c01a71f 340 fTOFTPCPtAxis->SetName("fTOFTPCPtAxis");
341 fTOFTPCPtAxis->SetTitle("p_{T} GeV/c");
07d62e30 342
343 // Create Pt spectrum histogram.
344 fPtSpectrumTOFTPCbins = new TH1F("fPtSpectrumTOFTPCbins","p_{T} Spectrum;p_{T} (GeV/c);Count",nptbins,ptarrayTOFTPC);
345 fOutputList->Add(fPtSpectrumTOFTPCbins);
346
347 // Create unidentified correlations histogram.
348 fCorrelationsTOFTPCbins = AliHistToolsDiHadronPID::MakeHist3D("fCorrelationsTOFTPCbins","Correlations;#Delta#phi;#Delta#eta;p_{T} (GeV/c)",
349 fNDPhiBins,-TMath::Pi()/2.,3.*TMath::Pi()/2.,
350 fNDEtaBins,-1.6,1.6,
351 nptbins, ptarrayTOFTPC);
352 fOutputList->Add(fCorrelationsTOFTPCbins);
353
354 // Create unidentified mixed events histogram.
355 fMixedEventsTOFTPCbins = AliHistToolsDiHadronPID::MakeHist3D("fMixedEventsTOFTPCbins","Mixed Events;#Delta#phi;#Delta#eta;p_{T} (GeV/c)",
356 fNDPhiBins,-TMath::Pi()/2.,3.*TMath::Pi()/2.,
357 fNDEtaBins,-1.6,1.6,
358 nptbins, ptarrayTOFTPC);
359 fOutputList->Add(fMixedEventsTOFTPCbins);
5c01a71f 360
5c01a71f 361 fTOFTPChistos = new TObjArray(3);
362 fTOFTPChistos->SetOwner(kTRUE);
07d62e30 363 fTOFTPChistos->SetName("CorrelationsTOFTPC");
364
365 if (fCalculateMismatch) {
366 fTOFTPCmismatch = new TObjArray(3);
367 fTOFTPCmismatch->SetOwner(kTRUE);
368 fTOFTPCmismatch->SetName("MismatchTOFTPC");
369 }
5c01a71f 370
371 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
372
0770be27 373 // Create the directory structure Pion, Kaon, Proton, regardless
374 // of wether the histograms are created (to keep the order.)
07d62e30 375 TObjArray* TOFTPChistosTmp = new TObjArray(fTOFTPCPtAxis->GetNbins());
376 TOFTPChistosTmp->SetOwner(kTRUE);
377 TOFTPChistosTmp->SetName(speciesname[iSpecies].Data());
378
379 TObjArray* TOFTPCmismatchTmp = 0x0;
380 if (fCalculateMismatch) {
381 TOFTPCmismatchTmp = new TObjArray(fTOFTPCPtAxis->GetNbins());
382 TOFTPCmismatchTmp->SetOwner(kTRUE);
383 TOFTPCmismatchTmp->SetName(speciesname[iSpecies].Data());
384 }
5c01a71f 385
0770be27 386 // Only Create the TOF/TPC histograms when requested.
387 Bool_t MakeTOFTPCcorrelations[3] = {fMakeTOFTPCcorrelationsPi, fMakeTOFTPCcorrelationsKa, fMakeTOFTPCcorrelationsPr};
388 if (MakeTOFTPCcorrelations[iSpecies]) {
389 for (Int_t iBinPt = 1; iBinPt < (fTOFTPCPtAxis->GetNbins() + 1); iBinPt++) {
390
391 // Approximate resolutions of TOF and TPC detector.
392 const Double_t sTOFest = 110.;
393 const Double_t sTPCest = 4.5;
394
395 // Set range +/- 5 sigma of main peak. (+ 10 sigma for TOF max, for mismatches.)
396 Double_t TOFmin = -5. * sTOFest;
397 Double_t TOFmax = 10. * sTOFest;
398 Double_t TPCmin = -4. * sTPCest;
399 Double_t TPCmax = 4. * sTPCest;
400
401 Double_t TOFexp = AliFunctionsDiHadronPID::TOFExpTime(fTOFTPCPtAxis->GetBinLowEdge(iBinPt), 0.4, AliFunctionsDiHadronPID::M(iSpecies));
402 Double_t TPCexp = AliFunctionsDiHadronPID::TPCExpdEdX(fTOFTPCPtAxis->GetBinLowEdge(iBinPt), 0.4, AliFunctionsDiHadronPID::M(iSpecies));
403
404 for (Int_t jSpecies = 0; jSpecies < 3; jSpecies++) {
405
406 if (iSpecies == jSpecies) {continue;}
407
408 Double_t TOFexpOther = AliFunctionsDiHadronPID::TOFExpTime(fTOFTPCPtAxis->GetBinLowEdge(iBinPt), 0.4, AliFunctionsDiHadronPID::M(jSpecies));
409 Double_t TPCexpOther = AliFunctionsDiHadronPID::TPCExpdEdX(fTOFTPCPtAxis->GetBinLowEdge(iBinPt), 0.4, AliFunctionsDiHadronPID::M(jSpecies));
410
411 // If any peak is within +/- 7 sigma, then also add this peak.
412 if ( (TMath::Abs(TOFexp - TOFexpOther) < 7. * sTOFest) ||
413 (TMath::Abs(TPCexp - TPCexpOther) < 7. * sTPCest) ) {
414
415 TOFmin = TMath::Min(TOFmin, (TOFexpOther - TOFexp - 2. * sTOFest) );
416 TOFmax = TMath::Max(TOFmax, (TOFexpOther - TOFexp + 10. * sTOFest) );
417 TPCmin = TMath::Min(TPCmin, (TPCexpOther - TPCexp - 2. * sTPCest) );
418 TPCmax = TMath::Max(TPCmax, (TPCexpOther - TPCexp + 2. * sTPCest) );
5c01a71f 419
0770be27 420 }
5c01a71f 421
422 }
423
0770be27 424 // With the standard TOF range, fitting the deuterons and the TOF mismatches is
425 // hard. This flag doubles the range of the TOF axis in the TOF/TPC histograms,
426 // while leaving the resolution the same. Turning on this flag will greatly increase
427 // the memory consumption of the task, to the point that it's probably too much
428 // to save a Buffer with all three species included.
429 Double_t TOFreach = TOFmax - TOFmin;
430 TOFmax += (TOFreach * (fTOFIntervalFactorTOFTPC - 1.));
431 Int_t TOFbins = (60. * fTOFIntervalFactorTOFTPC);
432
433 Int_t NBinsTOFTPC[4] = {32, 32, TOFbins, 40};
434 Double_t minTOFTPC[4] = {-TMath::Pi()/2., -1.6, TOFmin, TPCmin};
435 Double_t maxTOFTPC[4] = {3.*TMath::Pi()/2., 1.6, TOFmax, TPCmax};
436
437 THnF* htmp = new THnF(Form("fCorrelationsTOFTPC_%i",iBinPt),
438 Form("%5.3f < p_{T} < %5.3f", fTOFTPCPtAxis->GetBinLowEdge(iBinPt), fTOFTPCPtAxis->GetBinUpEdge(iBinPt)),
439 4, NBinsTOFTPC, minTOFTPC, maxTOFTPC);
440
441 (htmp->GetAxis(0))->SetTitle("#Delta#phi");
442 (htmp->GetAxis(1))->SetTitle("#Delta#eta");
443 (htmp->GetAxis(2))->SetTitle("t_{TOF} (ps)");
444 (htmp->GetAxis(3))->SetTitle("dE/dx (a.u.)");
445
446 TOFTPChistosTmp->Add(htmp);
447
448 if (fCalculateMismatch) {
449 // Mismatch histogram.
450 TH2F* htmp2 = new TH2F(Form("fMismatchTOFTPC_%i",iBinPt),
451 Form("%5.3f < p_{T} < %5.3f; t_{TOF} (ps); dE/dx (a.u.)", fTOFTPCPtAxis->GetBinLowEdge(iBinPt), fTOFTPCPtAxis->GetBinUpEdge(iBinPt)),
452 NBinsTOFTPC[2], TOFmin, TOFmax, NBinsTOFTPC[3], TPCmin, TPCmax);
453 htmp2->SetDirectory(0);
454
455 TOFTPCmismatchTmp->Add(htmp2);
456
457 }
458 } // End loop over pT bins.
459 } // End species if.
5c01a71f 460
07d62e30 461 fTOFTPChistos->Add(TOFTPChistosTmp);
462 if (fCalculateMismatch) {fTOFTPCmismatch->Add(TOFTPCmismatchTmp);}
5c01a71f 463
464 }
465
466 fOutputList->Add(fTOFTPChistos);
07d62e30 467 if (fCalculateMismatch) {fOutputList->Add(fTOFTPCmismatch);}
5c01a71f 468
469 }
6788af99 470
471 // Load external TOF histograms if flag is set.
07d62e30 472 if (fCalculateMismatch) {LoadExtMismatchHistos();}
6788af99 473
474 PostData(1,fOutputList);
475
476}
477
97724bd1 478// -----------------------------------------------------------------------
6788af99 479void AliAnalysisTaskDiHadronPID::LocalInit() {
480
481 //
5c01a71f 482 // Initialize on the this computer.
6788af99 483 //
484
50dfda71 485 if (fDebug > 0) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
5c01a71f 486
6788af99 487}
488
97724bd1 489// -----------------------------------------------------------------------
6788af99 490void AliAnalysisTaskDiHadronPID::UserExec(Option_t*) {
491
492 //
493 // Main Loop.
494 //
495
50dfda71 496 if (fDebug > 0) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
6788af99 497
498 // Input Current Event.
499 fCurrentAODEvent = dynamic_cast<AliAODEvent*>(InputEvent());
500 if (!fCurrentAODEvent) AliFatal("No Event Found!");
501
502 if (!fEventCuts->IsSelected(fCurrentAODEvent)) {return;}
503
504 // Fill the global tracks array. - NOT NEEDED I THINK, since we're not using
505 // bit 1<<7 for the associated tracks!
506
507 // Let the track cut objects know that a new event will start.
508 fTrackCutsTrigger->StartNewEvent();
509 fTrackCutsAssociated->StartNewEvent();
510
511 // Create arrays for trigger/associated particles.
512 fTriggerTracks = new TObjArray(350);
513 fTriggerTracks->SetOwner(kTRUE);
514
515 fAssociatedTracks = new TObjArray(3500);
516 fAssociatedTracks->SetOwner(kTRUE);
517
518 for (Int_t iTrack = 0; iTrack < fCurrentAODEvent->GetNumberOfTracks(); iTrack++) {
519
520 AliAODTrack* track = (AliAODTrack*)fCurrentAODEvent->GetTrack(iTrack);
521 AliTrackDiHadronPID* pidtrack = new AliTrackDiHadronPID(track,0x0,0x0,fPIDResponse);
522 pidtrack->ForgetAboutPointers();
523 pidtrack->SetDebugLevel(fDebug);
524
525 Double_t rndhittime = -1.e21;
07d62e30 526 if (fCalculateMismatch) rndhittime = GenerateRandomHit(pidtrack->Eta());
6788af99 527
528 // Fill the trigger/associated tracks array.
529 if (fTrackCutsTrigger->IsSelectedData(pidtrack,rndhittime)) {fTriggerTracks->AddLast(pidtrack);}
07d62e30 530 else if (fTrackCutsAssociated->IsSelectedData(pidtrack,rndhittime)) {
531
532 fAssociatedTracks->AddLast(pidtrack);
533
534 // Fill p_T spectrum.
535 if (fPtSpectrumTOFbins) fPtSpectrumTOFbins->Fill(pidtrack->Pt());
536 if (fPtSpectrumTOFTPCbins) fPtSpectrumTOFTPCbins->Fill(pidtrack->Pt());
537
538 // Fill mismatch histograms with associateds.
539 if (fCalculateMismatch && (rndhittime > -1.e20)) {
540
541 Double_t apt = pidtrack->Pt();
542
543 if (fMakeTOFcorrelations) {
544
545 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
546
547 TObjArray* atmp = (TObjArray*)fTOFmismatch->At(iSpecies);
548 Int_t ptbin = fTOFPtAxis->FindBin(apt);
549
550 // Only fill if histogram exists in fTOFmismatch.
551 if ( !(ptbin < 1) && !(ptbin > fTOFPtAxis->GetNbins()) ) {
552
553 TH1F* htmp = (TH1F*)atmp->At(ptbin - 1);
554 htmp->Fill(rndhittime - pidtrack->GetTOFsignalExpected(iSpecies));
555
556 }
557 }
558 }
559
0770be27 560 Bool_t MakeTOFTPCcorrelations[3] = {fMakeTOFTPCcorrelationsPi, fMakeTOFTPCcorrelationsKa, fMakeTOFTPCcorrelationsPr};
561 if (fMakeTOFTPCcorrelationsPi || fMakeTOFTPCcorrelationsKa || fMakeTOFTPCcorrelationsPr) {
07d62e30 562
563 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
564
0770be27 565 if (!MakeTOFTPCcorrelations[iSpecies]) {continue;}
566
07d62e30 567 TObjArray* atmp = (TObjArray*)fTOFTPCmismatch->At(iSpecies);
568 Int_t ptbin = fTOFTPCPtAxis->FindBin(apt);
569
570 // Only fill if histogram exists in fTOFTPCmismatch.
571 if ( !(ptbin < 1) && !(ptbin > fTOFTPCPtAxis->GetNbins()) ) {
572
573 TH2F* htmp = (TH2F*)atmp->At(ptbin - 1);
574 htmp->Fill(rndhittime - pidtrack->GetTOFsignalExpected(iSpecies), pidtrack->GetTPCsignalMinusExpected(iSpecies));
575
576 }
577 }
578 }
579
580 }
581
582 }
6788af99 583 else {delete pidtrack;}
584
585 }
586
587 // Fill Correlation histograms.
588 for (Int_t iTrigger = 0; iTrigger < fTriggerTracks->GetEntriesFast(); iTrigger++) {
589 AliTrackDiHadronPID* triggertrack = (AliTrackDiHadronPID*)fTriggerTracks->At(iTrigger);
590
591 for (Int_t iAssociated = 0; iAssociated < fAssociatedTracks->GetEntriesFast(); iAssociated++) {
592 AliTrackDiHadronPID* associatedtrack = (AliTrackDiHadronPID*)fAssociatedTracks->At(iAssociated);
593
594 Double_t DPhi = triggertrack->Phi() - associatedtrack->Phi();
595 if (DPhi < -TMath::Pi()/2.) {DPhi += 2.*TMath::Pi();}
596 if (DPhi > 3.*TMath::Pi()/2.) {DPhi -= 2.*TMath::Pi();}
597
598 Double_t DEta = triggertrack->Eta() - associatedtrack->Eta();
07d62e30 599 if (fCorrelationsTOFbins) fCorrelationsTOFbins->Fill(DPhi,DEta,associatedtrack->Pt());
600 if (fCorrelationsTOFTPCbins) fCorrelationsTOFTPCbins->Fill(DPhi,DEta,associatedtrack->Pt());
6788af99 601
0770be27 602 Double_t apt = associatedtrack->Pt();
5c01a71f 603
0770be27 604 // Fill TOF correlations.
605 if (fMakeTOFcorrelations) {
5c01a71f 606
0770be27 607 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
608
5c01a71f 609 TObjArray* atmp = (TObjArray*)fTOFhistos->At(iSpecies);
610 Int_t ptbin = fTOFPtAxis->FindBin(apt);
6788af99 611
5c01a71f 612 // Only fill if histogram exists in fTOFhistos.
613 if ( !(ptbin < 1) && !(ptbin > fTOFPtAxis->GetNbins()) ) {
614
615 TH3F* htmp = (TH3F*)atmp->At(ptbin - 1);
616 htmp->Fill(DPhi, DEta, associatedtrack->GetTOFsignalMinusExpected(iSpecies));
6788af99 617
6788af99 618 }
5c01a71f 619 }
0770be27 620 }
621
622 // Fill TOF/ TPC Correlations.
623 Bool_t MakeTOFTPCcorrelations[3] = {fMakeTOFTPCcorrelationsPi, fMakeTOFTPCcorrelationsKa, fMakeTOFTPCcorrelationsPr};
624 if (fMakeTOFTPCcorrelationsPi || fMakeTOFTPCcorrelationsKa || fMakeTOFTPCcorrelationsPr) {
5c01a71f 625
0770be27 626 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
627
628 if (!MakeTOFTPCcorrelations[iSpecies]) {continue;}
5c01a71f 629
630 TObjArray* atmp = (TObjArray*)fTOFTPChistos->At(iSpecies);
631 Int_t ptbin = fTOFTPCPtAxis->FindBin(apt);
632
633 // Only fill if histogram exists in fTOFhistos.
634 if ( !(ptbin < 1) && !(ptbin > fTOFTPCPtAxis->GetNbins()) ) {
635
636 THnF* htmp = (THnF*)atmp->At(ptbin - 1);
637 Double_t TOFTPCfill[4] = {DPhi, DEta,
638 associatedtrack->GetTOFsignalMinusExpected(iSpecies), associatedtrack->GetTPCsignalMinusExpected(iSpecies)};
639
640 htmp->Fill(TOFTPCfill);
6788af99 641
5c01a71f 642 }
0770be27 643 }
644 }
6788af99 645 }
646 }
647
07d62e30 648 //cout<<"Triggers: "<<fTriggerTracks->GetEntriesFast()<<" Associateds: "<<fAssociatedTracks->GetEntriesFast()<<endl;
6788af99 649
650 // Determine vtxz of current event.
651 AliAODVertex* currentprimaryvertex = fCurrentAODEvent->GetPrimaryVertex();
652 Double_t vtxz = currentprimaryvertex->GetZ();
653
07d62e30 654 // Determine centrality of current event (for PbPb).
655 AliEventPool* poolin = 0x0;
656 Float_t percentile = -1.;
657 if (fEventCuts->GetIsPbPb()) {
658 TString centralityestimator = fEventCuts->GetCentralityEstimator();
659 AliCentrality* currentcentrality = fCurrentAODEvent->GetCentrality();
660 percentile = currentcentrality->GetCentralityPercentile(centralityestimator.Data());
661
662 poolin = fPoolMgr->GetEventPool(percentile, vtxz);
663 if (!poolin) {AliFatal(Form("No pool found for centrality = %f, vtxz = %f", percentile, vtxz));}
664 } else {
665 poolin = fPoolMgr->GetEventPool(0.5, vtxz); // There are no multiplicity bins for pp yet.
666 if (!poolin) {AliFatal(Form("No pool found for vtxz = %f", vtxz));}
667 }
668
6788af99 669 // TObjArray* fGlobalTracksArray;
670
50dfda71 671 // Give a print out of the pool manager's contents.
f054df96 672 if (fDebug > 0) PrintPoolManagerContents();
50dfda71 673
674 // Mix events if there are enough events in the pool.
6788af99 675 if (poolin->GetCurrentNEvents() >= fMinNEventsForMixing) {
07d62e30 676 //{cout << "Mixing Events." << endl;}
50dfda71 677
678 // Loop over all events in the event pool.
679 for (Int_t iMixEvent = 0; iMixEvent < poolin->GetCurrentNEvents(); iMixEvent++) {
680 TObjArray* mixtracks = poolin->GetEvent(iMixEvent);
681
682 // Mix either the triggers or the associateds.
683 if (fMixTriggers) {
684
685 // Loop over all associateds in this event.
686 for (Int_t iAssociated = 0; iAssociated < fAssociatedTracks->GetEntriesFast(); iAssociated++) {
687 AliTrackDiHadronPID* associatedtrack = (AliTrackDiHadronPID*)fAssociatedTracks->At(iAssociated);
688
689 // Loop over all mixed tracks.
690 for (Int_t iMixTrack = 0; iMixTrack < mixtracks->GetEntriesFast(); iMixTrack++) {
691 AliTrackDiHadronPID* mixtrack = (AliTrackDiHadronPID*)mixtracks->At(iMixTrack);
692
693 Double_t DPhi = mixtrack->Phi() - associatedtrack->Phi();
694 if (DPhi < -TMath::Pi()/2.) {DPhi += 2.*TMath::Pi();}
695 if (DPhi > 3.*TMath::Pi()/2.) {DPhi -= 2.*TMath::Pi();}
696
697 Double_t DEta = mixtrack->Eta() - associatedtrack->Eta();
07d62e30 698 if (fMixedEventsTOFbins) fMixedEventsTOFbins->Fill(DPhi,DEta,associatedtrack->Pt());
699 if (fMixedEventsTOFTPCbins) fMixedEventsTOFTPCbins->Fill(DPhi,DEta,associatedtrack->Pt());
50dfda71 700
701 }
702 }
6788af99 703
50dfda71 704 } else {
6788af99 705
50dfda71 706 // Loop over all triggers in this event.
707 for (Int_t iTrigger = 0; iTrigger < fTriggerTracks->GetEntriesFast(); iTrigger++) {
708 AliTrackDiHadronPID* triggertrack = (AliTrackDiHadronPID*)fTriggerTracks->At(iTrigger);
6788af99 709
50dfda71 710 // Loop over all mixed tracks.
711 for (Int_t iMixTrack = 0; iMixTrack < mixtracks->GetEntriesFast(); iMixTrack++) {
712 AliTrackDiHadronPID* mixtrack = (AliTrackDiHadronPID*)mixtracks->At(iMixTrack);
713
714 Double_t DPhi = triggertrack->Phi() - mixtrack->Phi();
715 if (DPhi < -TMath::Pi()/2.) {DPhi += 2.*TMath::Pi();}
716 if (DPhi > 3.*TMath::Pi()/2.) {DPhi -= 2.*TMath::Pi();}
6788af99 717
50dfda71 718 Double_t DEta = triggertrack->Eta() - mixtrack->Eta();
07d62e30 719 if (fMixedEventsTOFbins) fMixedEventsTOFbins->Fill(DPhi,DEta,mixtrack->Pt());
720 if (fMixedEventsTOFTPCbins) fMixedEventsTOFTPCbins->Fill(DPhi,DEta,mixtrack->Pt());
50dfda71 721 }
722 }
723
724 } // End if
6788af99 725 }
726 }
727
728 // Update the event pool.
07d62e30 729 AliEventPool* poolout = 0x0;
730 if (fEventCuts->GetIsPbPb()) {
731 poolout = fPoolMgr->GetEventPool(percentile, vtxz); // Get the buffer associated with the current centrality and z-vtx
732 if (!poolout) AliFatal(Form("No pool found for centrality = %f, vtx_z = %f", percentile, vtxz));
733 } else {
734 poolout = fPoolMgr->GetEventPool(0.5, vtxz); // Get the buffer associated with the current centrality and z-vtx
735 if (!poolout) AliFatal(Form("No pool found for vtx_z = %f", vtxz));
736 }
737
6788af99 738
739 // Q: is it a problem that the fAssociatedTracks array can be bigger than the number of tracks inside?
50dfda71 740 if (fMixTriggers) {
741 poolout->UpdatePool(fTriggerTracks);
742 fAssociatedTracks->Delete();
743 delete fAssociatedTracks;
744 }
745 else {
746 poolout->UpdatePool(fAssociatedTracks);
747 fTriggerTracks->Delete();
748 delete fTriggerTracks;
749 }
6788af99 750
6788af99 751 fTriggerTracks = 0x0;
752 fAssociatedTracks = 0x0;
753
754 // Tell the track cut object that the event is done.
755 fTrackCutsTrigger->EventIsDone(0);
756 fTrackCutsAssociated->EventIsDone(0);
757
758 PostData(1,fOutputList);
759
760}
761
97724bd1 762// -----------------------------------------------------------------------
6788af99 763Bool_t AliAnalysisTaskDiHadronPID::LoadExtMismatchHistos() {
764
765 //
766 // Attempting to load a root file containing information needed
767 // to generate random TOF hits.
768 //
769
50dfda71 770 if (fDebug > 0) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
6788af99 771
772 // Opening external TOF file.
773 if (fDebug > 0) cout<<"Trying to open TOFmismatchHistos.root ..."<<endl;
774 TFile* fin = 0x0;
775 fin = TFile::Open("alien:///alice/cern.ch/user/m/mveldhoe/rootfiles/TOFmismatchHistos.root");
776 if (!fin) {
777 AliWarning("Couln't open TOFmismatchHistos, will not calculate mismatches...");
07d62e30 778 fCalculateMismatch = kFALSE;
6788af99 779 return kFALSE;
780 }
781
782 // Check if the required histograms are present.
783 TH1F* tmp1 = (TH1F*)fin->Get("hNewT0Fill");
784 if (!tmp1) {
785 AliWarning("Couln't find hNewT0Fill, will not calculate mismatches...");
07d62e30 786 fCalculateMismatch = kFALSE;
6788af99 787 return kFALSE;
788 }
789 TH2F* tmp2 = (TH2F*)fin->Get("hLvsEta");
790 if (!tmp2) {
791 AliWarning("Couln't find hLvsEta, will not calculate mismatches...");
07d62e30 792 fCalculateMismatch = kFALSE;
6788af99 793 return kFALSE;
794 }
795
796 // Make a deep copy of the files in the histogram.
797 fT0Fill = (TH1F*)tmp1->Clone("fT0Fill");
798 fLvsEta = (TH2F*)tmp2->Clone("fLvsEta");
799
800 // Close the external file.
801 AliInfo("Closing external file.");
802 fin->Close();
803
804 // Creating a TObjArray for LvsEta projections.
805 const Int_t nbinseta = fLvsEta->GetNbinsX();
806 fLvsEtaProjections = new TObjArray(nbinseta);
807 fLvsEtaProjections->SetOwner(kTRUE);
808
809 // Making the projections needed (excluding underflow/ overflow).
810 for (Int_t iEtaBin = 1; iEtaBin < (nbinseta + 1); iEtaBin++) {
811 TH1F* tmp = (TH1F*)fLvsEta->ProjectionY(Form("LvsEtaProjection_%i",iEtaBin),iEtaBin,iEtaBin);
812 tmp->SetDirectory(0);
07d62e30 813 fLvsEtaProjections->AddAt(tmp,iEtaBin - 1);
6788af99 814 }
815
816 return kTRUE;
817
818}
819
97724bd1 820// -----------------------------------------------------------------------
6788af99 821Double_t AliAnalysisTaskDiHadronPID::GenerateRandomHit(Double_t eta) {
822
823 //
824 // Returns a random TOF time.
825 //
826
50dfda71 827 if (fDebug > 0) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
6788af99 828
829 // Default (error) value:
830 Double_t rndhittime = -1.e21;
831
832 // TOF mismatch flag is not turned on.
07d62e30 833 if (!fCalculateMismatch) {
834 AliFatal("Called GenerateRandomHit() method, but flag fCalculateMismatch not set.");
6788af99 835 return rndhittime;
836 }
837
838 // TOF doesn't extend much further than 0.8.
839 if (TMath::Abs(eta) > 0.8) {
840 if (fDebug) {AliInfo("Tried to get a random hit for a track with eta > 0.8.");}
841 return rndhittime;
842 }
843
844 // Finding the bin of the eta.
845 TAxis* etaAxis = fLvsEta->GetXaxis();
846 Int_t etaBin = etaAxis->FindBin(eta);
07d62e30 847 if (etaBin == 0 || (etaBin == etaAxis->GetNbins() + 1)) {return rndhittime;}
848
849 const TH1F* lengthDistribution = (const TH1F*)fLvsEtaProjections->At(etaBin - 1);
6788af99 850
851 if (!lengthDistribution) {
852 AliFatal("length Distribution not found.");
853 return rndhittime;
854 }
855
856 Double_t currentRndLength = lengthDistribution->GetRandom(); // in cm.
857
858 // Similar to Roberto's code.
859 Double_t currentRndTime = currentRndLength / (TMath::C() * 1.e2 / 1.e12);
860 Double_t t0fill = -1.26416e+04;
861 rndhittime = fT0Fill->GetRandom() - t0fill + currentRndTime;
862
863 return rndhittime;
864
865}
866
97724bd1 867// -----------------------------------------------------------------------
50dfda71 868void AliAnalysisTaskDiHadronPID::PrintPoolManagerContents() {
869
870 //
871 // Prints out the current contents of the event pool manager.
872 //
873
874 // Determine the number of pools in the pool manager.
875 AliEventPool* poolin = fPoolMgr->GetEventPool(0,0);
876 Int_t NPoolsCentrality = 0;
877 while (poolin) {
878 NPoolsCentrality++;
879 poolin = fPoolMgr->GetEventPool(NPoolsCentrality,0);
880 }
881
882 poolin = fPoolMgr->GetEventPool(0,0);
883 Int_t NPoolsVtxZ = 0;
884 while (poolin) {
885 NPoolsVtxZ++;
886 poolin = fPoolMgr->GetEventPool(0,NPoolsVtxZ);
887 }
888
889 // Loop over all Pools in the matrix of the pool manager.
890 cout<<" Pool manager contents: (Nevt,NTrack)"<<endl;
891 for (Int_t iCentrality = 0; iCentrality < NPoolsCentrality; iCentrality++) {
892 cout<<Form("Centrality Bin: %2i --> ", iCentrality);
893
894 for (Int_t iVtxZ = 0; iVtxZ < NPoolsVtxZ; iVtxZ++) {
895
896 poolin = fPoolMgr->GetEventPool(iCentrality, iVtxZ);
897
898 cout<<Form("(%2i,%4i) ",poolin->GetCurrentNEvents(), poolin->NTracksInPool());
899
900 }
901
902 cout<<endl;
903 }
904
905}
906
97724bd1 907// -----------------------------------------------------------------------
6788af99 908void AliAnalysisTaskDiHadronPID::Terminate(Option_t*) {;
909
910 //
911 // Called when task is done.
912 //
913
50dfda71 914 if (fDebug > 0) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
6788af99 915
916 delete fT0Fill;
917 fT0Fill = 0x0;
918 delete fLvsEta;
919 fLvsEta = 0x0;
920 delete fLvsEtaProjections;
921 fLvsEtaProjections = 0x0;
922
923}