]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGJE/UserTasks/AliAnalysisTaskJetProtonCorr.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
add common abstract interface classes for flat and fat ESDs
[u/mrichter/AliRoot.git] / PWGJE / UserTasks / AliAnalysisTaskJetProtonCorr.cxx
... / ...
CommitLineData
1// ROOT
2#include "TFile.h"
3#include "TList.h"
4#include "TH1.h"
5#include "TH2.h"
6#include "TH3.h"
7#include "TF1.h"
8#include "TFormula.h"
9#include "TRandom.h"
10
11// analysis framework
12#include "AliAnalysisManager.h"
13#include "AliInputEventHandler.h"
14#include "AliVEvent.h"
15#include "AliVTrack.h"
16#include "AliVTrdTrack.h"
17#include "AliVVertex.h"
18#include "AliPIDResponse.h"
19#include "AliEventPoolManager.h"
20#include "AliOADBContainer.h"
21#include "AliTOFPIDParams.h"
22#include "AliAnalysisTaskVnV0.h"
23
24// MC stuff
25#include "AliMCEvent.h"
26#include "AliGenPythiaEventHeader.h"
27
28// ESD stuff
29#include "AliESDEvent.h"
30#include "AliESDInputHandler.h"
31#include "AliESDtrack.h"
32#include "AliESDtrackCuts.h"
33#include "AliESDTrdTrack.h"
34#include "AliESDTrdTracklet.h"
35#include "AliESDTrdTrigger.h"
36
37// AOD stuff
38#include "AliAODEvent.h"
39#include "AliAODJet.h"
40#include "AliAODTrack.h"
41
42// jet tasks
43#include "AliAnalysisTaskJetServices.h"
44#include "AliAnalysisHelperJetTasks.h"
45
46#include "AliAnalysisTaskJetProtonCorr.h"
47
48#include <iostream>
49#include <cmath>
50
51AliAnalysisTaskJetProtonCorr::AliAnalysisTaskJetProtonCorr(const char *name) :
52 AliAnalysisTaskSE(name),
53 fMCEvent(0x0),
54 fESDEvent(0x0),
55 fAODEvent(0x0),
56 fRunNumber(-1),
57 fOADBContainerTOF(0x0),
58 fParamsTOF(0x0),
59 fEventplane(0x0),
60 fTriggerMask(0),
61 fClassMask(0),
62 fCentrality(100.),
63 fCentralityCheck(100.),
64 fZvtx(0.),
65 fPIDResponse(0x0),
66 fEventPlaneAngle(5.),
67 fEventPlaneAngleCheck(5.),
68 fEventPlaneAngle3(5.),
69 fPrimTrackArrayAss(0x0),
70 fPrimTrackArrayTrg(0x0),
71 fPrimConstrainedTrackArray(new TClonesArray("AliESDtrack", 100)),
72 fJetArray(0x0),
73 fPoolMgr(),
74 fPool(),
75 fHistCorr(0x0),
76 fErrorMsg(10),
77 fOutputList(),
78 fHist(),
79 fShortTaskId("jet_prot_corr"),
80 fUseStandardCuts(kTRUE),
81 fUseEvplaneV0(kFALSE),
82 fCutsPrimTrg(0x0),
83 fCutsPrimTrgConstrain(new AliESDtrackCuts()),
84 fCutsPrimAss(0x0),
85 fCutsTwoTrackEff(0.02),
86 fAssFilterMask(1 << 10),
87 fRequirePID(kTRUE),
88 fTrgJetEtaMax(0.45),
89 fHadEtaMax(0.8),
90 fTrgPartPtMin(6.),
91 fTrgPartPtMax(8.),
92 fTrgJetPtMin(50.),
93 fTrgJetPtMax(80.),
94 fTrgJetLeadTrkPtMin(6.),
95 fTrgJetLeadTrkPtMax(100.),
96 fTrgJetAreaMin(0.6 * TMath::Pi() * 0.2*0.2),
97 fAssPartPtMin(2.),
98 fAssPartPtMax(4.),
99 fTrgAngleToEvPlane(TMath::Pi() / 4.),
100 fToyMeanNoPart(1.),
101 fToyRadius(.8),
102 fToySmearPhi(.2),
103 fTrgJetPhiModCent(new TF1("jetphimodcent", "1 + 2 * [0] * cos(2*x)", 0., 2 * TMath::Pi())),
104 fTrgJetPhiModSemi(new TF1("jetphimodsemi", "1 + 2 * [0] * cos(2*x)", 0., 2 * TMath::Pi())),
105 fTrgHadPhiModCent(new TF1("hadphimodcent", "1 + 2 * [0] * cos(2*x)", 0., 2 * TMath::Pi())),
106 fTrgHadPhiModSemi(new TF1("hadphimodsemi", "1 + 2 * [0] * cos(2*x)", 0., 2 * TMath::Pi()))
107{
108 // default ctor
109
110 fkCorrTypeName[kCorrHadHad] = "hh";
111 fkCorrTypeName[kCorrHadProt] = "hp";
112 fkCorrTypeName[kCorrJetHad] = "jh";
113 fkCorrTypeName[kCorrJetProt] = "jp";
114 fkCorrTypeName[kCorrRndJetHad] = "rjh";
115 fkCorrTypeName[kCorrRndJetProt] = "rjp";
116 fkCorrTypeName[kCorrRndHadHad] = "rhh";
117 fkCorrTypeName[kCorrRndHadProt] = "rhp";
118 fkCorrTypeName[kCorrRndJetExcHad] = "rjeh";
119 fkCorrTypeName[kCorrRndJetExcProt] = "rjep";
120 fkCorrTypeName[kCorrRndHadExcHad] = "rheh";
121 fkCorrTypeName[kCorrRndHadExcProt] = "rhep";
122
123 fkClassName[kClCentral] = "cent";
124 fkClassName[kClSemiCentral] = "semi";
125 // fkClassName[kClDijet] = "dijet";
126
127 fkEvName[kEvSame] = "same";
128 fkEvName[kEvMix] = "mixed";
129
130 // track cuts for associates
131 if (fUseStandardCuts) {
132 fCutsPrimAss = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
133 } else {
134 fCutsPrimAss = new AliESDtrackCuts();
135
136 // this is taken from PWGJE track cuts
137 TFormula *f1NClustersTPCLinearPtDep = new TFormula("f1NClustersTPCLinearPtDep","70.+30./20.*x");
138 fCutsPrimAss->SetMinNClustersTPCPtDep(f1NClustersTPCLinearPtDep,20.);
139 fCutsPrimAss->SetMinNClustersTPC(70);
140 fCutsPrimAss->SetMaxChi2PerClusterTPC(4);
141 fCutsPrimAss->SetRequireTPCStandAlone(kTRUE); //cut on NClustersTPC and chi2TPC Iter1
142 fCutsPrimAss->SetAcceptKinkDaughters(kFALSE);
143 fCutsPrimAss->SetRequireTPCRefit(kTRUE);
144 fCutsPrimAss->SetMaxFractionSharedTPCClusters(0.4);
145 // ITS
146 fCutsPrimAss->SetRequireITSRefit(kTRUE);
147 //accept secondaries
148 fCutsPrimAss->SetMaxDCAToVertexXY(2.4);
149 fCutsPrimAss->SetMaxDCAToVertexZ(3.2);
150 fCutsPrimAss->SetDCAToVertex2D(kTRUE);
151 //reject fakes
152 fCutsPrimAss->SetMaxChi2PerClusterITS(36);
153 fCutsPrimAss->SetMaxChi2TPCConstrainedGlobal(36);
154
155 fCutsPrimAss->SetRequireSigmaToVertex(kFALSE);
156
157 fCutsPrimAss->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
158 }
159
160 fCutsPrimAss->SetEtaRange(-fHadEtaMax, fHadEtaMax);
161 fCutsPrimAss->SetPtRange(0.15, 1E+15);
162
163 // track cuts for triggers
164 fCutsPrimTrg = new AliESDtrackCuts(*fCutsPrimAss);
165
166 // azimuthal modulation for triggers
167 fTrgJetPhiModCent->SetParameter(0, .10);
168 fTrgJetPhiModSemi->SetParameter(0, .10);
169 fTrgHadPhiModCent->SetParameter(0, .04);
170 fTrgHadPhiModSemi->SetParameter(0, .10);
171
172 // event mixing pool
173 Double_t centralityBins[] = {
174 0., 2., 4., 6., 8., 10., // central
175 30., 32., 34., 36., 38., 40., 42., 44., 46., 48., 50., // semi-central
176 90.
177 };
178 Int_t nCentralityBins = sizeof(centralityBins)/sizeof(centralityBins[0]);
179
180 Double_t vertexBins[] = {
181 -10., -8., -6., -4., -2., 0., 2., 4., 6., 8., 10.
182 };
183 Int_t nVertexBins = sizeof(vertexBins)/sizeof(vertexBins[0]);
184
185 Double_t psiBins[12];
186 Int_t nPsiBins = sizeof(psiBins)/sizeof(psiBins[0]);
187 for (Int_t iBin = 0; iBin < nPsiBins; ++iBin)
188 psiBins[iBin] = iBin * TMath::Pi()/nPsiBins;
189
190 const Int_t poolSize = 10; // unused by current event pool implementation
191 const Int_t trackDepth = 10000; // number of tracks to maintain in mixing buffer
192 const Float_t trackDepthFraction = 0.1;
193 const Int_t targetEvents = 1;
194 for (Int_t iAss = 0; iAss <= kAssProt; ++iAss) {
195 GetPoolMgr((Ass_t) iAss) =
196 new AliEventPoolManager(poolSize, trackDepth,
197 nCentralityBins, centralityBins,
198 nVertexBins, vertexBins);
199 GetPoolMgr((Ass_t) iAss)->SetTargetValues(trackDepth, trackDepthFraction, targetEvents);
200 }
201
202 fHistCorr = new AliHistCorr*[kEvLast*kCorrLast*kClLast];
203
204 DefineOutput(1, TList::Class());
205}
206
207AliAnalysisTaskJetProtonCorr::~AliAnalysisTaskJetProtonCorr()
208{
209 // dtor
210
211 // delete [] fHistCorr;
212}
213
214void AliAnalysisTaskJetProtonCorr::UserCreateOutputObjects()
215{
216 // create user output objects
217
218 // open OADB file and retrieve the OADB container with TOF parameters
219 TString oadbFileName =
220 TString::Format("%s/COMMON/PID/data/TOFPIDParams.root", AliAnalysisManager::GetOADBPath());
221 TFile *oadbFile = TFile::Open(oadbFileName);
222 if(!oadbFile->IsOpen())
223 AliFatal(Form("Cannot open OADB file %s", oadbFileName.Data()));
224 AliOADBContainer *oadbContainer =
225 (AliOADBContainer*) oadbFile->Get("TOFoadb");
226 if (!oadbContainer)
227 AliFatal("Cannot fetch OADB container for VZERO EP selection");
228 fOADBContainerTOF = new AliOADBContainer(*oadbContainer);
229 oadbFile->Close();
230 delete oadbFile;
231
232 // setup list
233 OpenFile(1);
234 fOutputList = new TList();
235 fOutputList->SetOwner();
236
237 // setup histograms
238 TH1 *hist;
239 TH1 *histStat = AddHistogram(ID(kHistStat), "event statistics;;counts",
240 kStatLast-1, .5, kStatLast-.5);
241 histStat->GetXaxis()->SetBinLabel(ID(kStatSeen));
242 histStat->GetXaxis()->SetBinLabel(ID(kStatTrg));
243 histStat->GetXaxis()->SetBinLabel(ID(kStatVtx));
244 histStat->GetXaxis()->SetBinLabel(ID(kStatEvCuts));
245 histStat->GetXaxis()->SetBinLabel(ID(kStatUsed));
246 histStat->GetXaxis()->SetBinLabel(ID(kStatCent));
247 histStat->GetXaxis()->SetBinLabel(ID(kStatEvPlane));
248 histStat->GetXaxis()->SetBinLabel(ID(kStatPID));
249 histStat->GetXaxis()->SetBinLabel(ID(kStatCentral));
250 histStat->GetXaxis()->SetBinLabel(ID(kStatSemiCentral));
251
252 AddHistogram(ID(kHistVertexNctb), "number of vertex contributors;N_{ctb};counts",
253 100, 0., 2000.);
254 AddHistogram(ID(kHistVertexZ), "z-position of primary vertex;z (cm);counts",
255 100, -50., 50.);
256
257 AddHistogram(ID(kHistCentrality), "centrality;C;counts",
258 110, -5., 105.);
259 hist = AddHistogram(ID(kHistCentralityUsed), "centrality used;C;event class",
260 110, -5., 105.,
261 kClLast, -.5, kClLast-.5);
262 hist->GetYaxis()->SetBinLabel(LAB(kClCentral));
263 hist->GetYaxis()->SetBinLabel(LAB(kClSemiCentral));
264 // hist->GetYaxis()->SetBinLabel(LAB(kClDijet));
265 AddHistogram(ID(kHistCentralityCheck), "centrality check;C;counts",
266 110, -5., 105.);
267 hist = AddHistogram(ID(kHistCentralityCheckUsed), "centrality check used;C;event class",
268 110, -5., 105.,
269 kClLast, -.5, kClLast-.5);
270 hist->GetYaxis()->SetBinLabel(LAB(kClCentral));
271 hist->GetYaxis()->SetBinLabel(LAB(kClSemiCentral));
272 // hist->GetYaxis()->SetBinLabel(LAB(kClDijet));
273 AddHistogram(ID(kHistCentralityVsMult), "centrality - multiplicity;centrality percentile (%);N_{prim}",
274 100, 0., 100.,
275 100, 0., 2500.);
276
277 AddHistogram(ID(kHistSignalTPC), "TPC dE/dx;p (GeV/#it{c});dE/dx (arb. units)",
278 100, 0., 10., 200, 0., 300.);
279 AddHistogram(ID(kHistSignalTOF), "TOF time of flight;p_{T} (GeV/#it{c});t (ns)",
280 100, 0., 10., 200, 0., 50.);
281 AddHistogram(ID(kHistBetaTOF), "TOF beta;p (GeV/#it{c}); #beta",
282 100, 0., 10.,
283 100, 0., 1.);
284
285 AddHistogram(ID(kHistDeltaTPC), "TPC dE/dx;p (GeV/#it{c});dE/dx (arb. units)",
286 100, 0., 10., 200, -100., 100.);
287 AddHistogram(ID(kHistDeltaTPCSemi), "TPC dE/dx;p (GeV/#it{c});dE/dx (arb. units)",
288 100, 0., 10., 200, -100., 100.);
289 AddHistogram(ID(kHistDeltaTOF), "TOF time of flight;p (GeV/#it{c});t (ns)",
290 100, 0., 10., 200, -2., 2.);
291 AddHistogram(ID(kHistDeltaTOFSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
292 100, 0., 10., 200, -2., 2.);
293
294 // Nsigma templates - central
295 AddHistogram(ID(kHistNsigmaTPCe), "TPC N_{#sigma,p} - e hypothesis;p (GeV/#it{c});N_{#sigma,p}",
296 100, 0., 10.,
297 100, -25., 25.);
298 AddHistogram(ID(kHistNsigmaTPCmu), "TPC N_{#sigma,p} - #mu hypothesis;p (GeV/#it{c});N_{#sigma,p}",
299 100, 0., 10.,
300 100, -25., 25.);
301 AddHistogram(ID(kHistNsigmaTPCpi), "TPC N_{#sigma,p} - #pi hypothesis;p (GeV/#it{c});N_{#sigma,p}",
302 100, 0., 10.,
303 100, -25., 25.);
304 AddHistogram(ID(kHistNsigmaTPCk), "TPC N_{#sigma,p} - K hypothesis;p (GeV/#it{c});N_{#sigma,p}",
305 100, 0., 10.,
306 100, -25., 25.);
307 AddHistogram(ID(kHistNsigmaTPCp), "TPC N_{#sigma,p} - p hypothesis;p (GeV/#it{c});N_{#sigma,p}",
308 100, 0., 10.,
309 100, -25., 25.);
310 AddHistogram(ID(kHistNsigmaTPCd), "TPC N_{#sigma,p} - d hypothesis;p (GeV/#it{c});N_{#sigma,p}",
311 100, 0., 10.,
312 100, -25., 25.);
313 AddHistogram(ID(kHistNsigmaTPCe_e), "TPC N_{#sigma,p} - e hypothesis (id. e);p (GeV/#it{c});N_{#sigma,p}",
314 100, 0., 10.,
315 100, -25., 25.);
316
317 AddHistogram(ID(kHistNsigmaTOFe), "TOF N_{#sigma,p} - e hypothesis;p (GeV/#it{c});N_{#sigma,p}",
318 100, 0., 10.,
319 200, -100., 100.);
320 AddHistogram(ID(kHistNsigmaTOFmu), "TOF N_{#sigma,p} - #mu hypothesis;p (GeV/#it{c});N_{#sigma,p}",
321 100, 0., 10.,
322 200, -100., 100.);
323 AddHistogram(ID(kHistNsigmaTOFpi), "TOF N_{#sigma,p} - #pi hypothesis;p (GeV/#it{c});N_{#sigma,p}",
324 100, 0., 10.,
325 200, -100., 100.);
326 AddHistogram(ID(kHistNsigmaTOFk), "TOF N_{#sigma,p} - K hypothesis;p (GeV/#it{c});N_{#sigma,p}",
327 100, 0., 10.,
328 200, -100., 100.);
329 AddHistogram(ID(kHistNsigmaTOFp), "TOF N_{#sigma,p} - p hypothesis;p (GeV/#it{c});N_{#sigma,p}",
330 100, 0., 10.,
331 200, -100., 100.);
332 AddHistogram(ID(kHistNsigmaTOFd), "TOF N_{#sigma,p} - d hypothesis;p (GeV/#it{c});N_{#sigma,p}",
333 100, 0., 10.,
334 200, -100., 100.);
335 AddHistogram(ID(kHistNsigmaTOFmismatch), "TOF N_{#sigma,p} - mismatch;p (GeV/#it{c});N_{#sigma,p}",
336 100, 0., 10.,
337 200, -100., 100.);
338
339 // Nsigma templates - semi-central
340 AddHistogram(ID(kHistNsigmaTPCeSemi), "TPC N_{#sigma,p} - e hypothesis;p (GeV/#it{c});N_{#sigma,p}",
341 100, 0., 10.,
342 100, -25., 25.);
343 AddHistogram(ID(kHistNsigmaTPCmuSemi), "TPC N_{#sigma,p} - #mu hypothesis;p (GeV/#it{c});N_{#sigma,p}",
344 100, 0., 10.,
345 100, -25., 25.);
346 AddHistogram(ID(kHistNsigmaTPCpiSemi), "TPC N_{#sigma,p} - #pi hypothesis;p (GeV/#it{c});N_{#sigma,p}",
347 100, 0., 10.,
348 100, -25., 25.);
349 AddHistogram(ID(kHistNsigmaTPCkSemi), "TPC N_{#sigma,p} - K hypothesis;p (GeV/#it{c});N_{#sigma,p}",
350 100, 0., 10.,
351 100, -25., 25.);
352 AddHistogram(ID(kHistNsigmaTPCpSemi), "TPC N_{#sigma,p} - p hypothesis;p (GeV/#it{c});N_{#sigma,p}",
353 100, 0., 10.,
354 100, -25., 25.);
355 AddHistogram(ID(kHistNsigmaTPCdSemi), "TPC N_{#sigma,p} - d hypothesis;p (GeV/#it{c});N_{#sigma,p}",
356 100, 0., 10.,
357 100, -25., 25.);
358 AddHistogram(ID(kHistNsigmaTPCe_eSemi), "TPC N_{#sigma,p} - e hypothesis (id. e);p (GeV/#it{c});N_{#sigma,p}",
359 100, 0., 10.,
360 100, -25., 25.);
361
362 AddHistogram(ID(kHistNsigmaTOFeSemi), "TOF N_{#sigma,p} - e hypothesis;p (GeV/#it{c});N_{#sigma,p}",
363 100, 0., 10.,
364 200, -100., 100.);
365 AddHistogram(ID(kHistNsigmaTOFmuSemi), "TOF N_{#sigma,p} - #mu hypothesis;p (GeV/#it{c});N_{#sigma,p}",
366 100, 0., 10.,
367 200, -100., 100.);
368 AddHistogram(ID(kHistNsigmaTOFpiSemi), "TOF N_{#sigma,p} - #pi hypothesis;p (GeV/#it{c});N_{#sigma,p}",
369 100, 0., 10.,
370 200, -100., 100.);
371 AddHistogram(ID(kHistNsigmaTOFkSemi), "TOF N_{#sigma,p} - K hypothesis;p (GeV/#it{c});N_{#sigma,p}",
372 100, 0., 10.,
373 200, -100., 100.);
374 AddHistogram(ID(kHistNsigmaTOFpSemi), "TOF N_{#sigma,p} - p hypothesis;p (GeV/#it{c});N_{#sigma,p}",
375 100, 0., 10.,
376 200, -100., 100.);
377 AddHistogram(ID(kHistNsigmaTOFdSemi), "TOF N_{#sigma,p} - d hypothesis;p (GeV/#it{c});N_{#sigma,p}",
378 100, 0., 10.,
379 200, -100., 100.);
380 AddHistogram(ID(kHistNsigmaTOFmismatchSemi), "TOF N_{#sigma,p} - mismatch;p (GeV/#it{c});N_{#sigma,p}",
381 100, 0., 10.,
382 200, -100., 100.);
383
384 // delta templates
385 AddHistogram(ID(kHistDeltaTOFe), "TOF #Delta;p (GeV/#it{c});t (ns)",
386 100, 0., 10., 200, -2., 2.);
387 AddHistogram(ID(kHistDeltaTOFmu), "TOF #Delta;p (GeV/#it{c});t (ns)",
388 100, 0., 10., 200, -2., 2.);
389 AddHistogram(ID(kHistDeltaTOFpi), "TOF #Delta;p (GeV/#it{c});t (ns)",
390 100, 0., 10., 200, -2., 2.);
391 AddHistogram(ID(kHistDeltaTOFk), "TOF #Delta;p (GeV/#it{c});t (ns)",
392 100, 0., 10., 200, -2., 2.);
393 AddHistogram(ID(kHistDeltaTOFp), "TOF #Delta;p (GeV/#it{c});t (ns)",
394 100, 0., 10., 200, -2., 2.);
395 AddHistogram(ID(kHistDeltaTOFd), "TOF #Delta;p (GeV/#it{c});t (ns)",
396 100, 0., 10., 200, -2., 2.);
397
398 AddHistogram(ID(kHistDeltaTOFeSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
399 100, 0., 10., 200, -2., 2.);
400 AddHistogram(ID(kHistDeltaTOFmuSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
401 100, 0., 10., 200, -2., 2.);
402 AddHistogram(ID(kHistDeltaTOFpiSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
403 100, 0., 10., 200, -2., 2.);
404 AddHistogram(ID(kHistDeltaTOFkSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
405 100, 0., 10., 200, -2., 2.);
406 AddHistogram(ID(kHistDeltaTOFpSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
407 100, 0., 10., 200, -2., 2.);
408 AddHistogram(ID(kHistDeltaTOFdSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
409 100, 0., 10., 200, -2., 2.);
410
411 // sigma comparisons
412 // AddHistogram(ID(kHistExpSigmaTOFe), "TOF time of flight;p (GeV/#it{c});t (ns)",
413 // 100, 0., 10., 200, 0., .25);
414 // AddHistogram(ID(kHistExpSigmaTOFmu), "TOF time of flight;p (GeV/#it{c});t (ns)",
415 // 100, 0., 10., 200, 0., .25);
416 // AddHistogram(ID(kHistExpSigmaTOFpi), "TOF time of flight;p (GeV/#it{c});t (ns)",
417 // 100, 0., 10., 200, 0., .25);
418 // AddHistogram(ID(kHistExpSigmaTOFk), "TOF time of flight;p (GeV/#it{c});t (ns)",
419 // 100, 0., 10., 200, 0., .25);
420 // AddHistogram(ID(kHistExpSigmaTOFp), "TOF time of flight;p (GeV/#it{c});t (ns)",
421 // 100, 0., 10., 200, 0., .25);
422 // AddHistogram(ID(kHistExpSigmaTOFd), "TOF time of flight;p (GeV/#it{c});t (ns)",
423 // 100, 0., 10., 200, 0., .25);
424
425 // AddHistogram(ID(kHistExpSigmaTOFeSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
426 // 100, 0., 10., 200, 0., .25);
427 // AddHistogram(ID(kHistExpSigmaTOFmuSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
428 // 100, 0., 10., 200, 0., .25);
429 // AddHistogram(ID(kHistExpSigmaTOFpiSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
430 // 100, 0., 10., 200, 0., .25);
431 // AddHistogram(ID(kHistExpSigmaTOFkSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
432 // 100, 0., 10., 200, 0., .25);
433 // AddHistogram(ID(kHistExpSigmaTOFpSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
434 // 100, 0., 10., 200, 0., .25);
435 // AddHistogram(ID(kHistExpSigmaTOFdSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
436 // 100, 0., 10., 200, 0., .25);
437
438 // AddHistogram(ID(kHistCmpSigmaTOFe), "#sigma comparison;exp #sigma;template #sigma",
439 // 200, 0., .25, 200, 0., .25);
440 // AddHistogram(ID(kHistCmpSigmaTOFmu), "#sigma comparison;exp #sigma;template #sigma",
441 // 200, 0., .25, 200, 0., .25);
442 // AddHistogram(ID(kHistCmpSigmaTOFpi), "#sigma comparison;exp #sigma;template #sigma",
443 // 200, 0., .25, 200, 0., .25);
444 // AddHistogram(ID(kHistCmpSigmaTOFk), "#sigma comparison;exp #sigma;template #sigma",
445 // 200, 0., .25, 200, 0., .25);
446 // AddHistogram(ID(kHistCmpSigmaTOFp), "#sigma comparison;exp #sigma;template #sigma",
447 // 200, 0., .25, 200, 0., .25);
448 // AddHistogram(ID(kHistCmpSigmaTOFd), "#sigma comparison;exp #sigma;template #sigma",
449 // 200, 0., .25, 200, 0., .25);
450
451 // AddHistogram(ID(kHistCmpSigmaTOFeSemi), "#sigma comparison;exp #sigma;template #sigma",
452 // 200, 0., .25, 200, 0., .25);
453 // AddHistogram(ID(kHistCmpSigmaTOFmuSemi), "#sigma comparison;exp #sigma;template #sigma",
454 // 200, 0., .25, 200, 0., .25);
455 // AddHistogram(ID(kHistCmpSigmaTOFpiSemi), "#sigma comparison;exp #sigma;template #sigma",
456 // 200, 0., .25, 200, 0., .25);
457 // AddHistogram(ID(kHistCmpSigmaTOFkSemi), "#sigma comparison;exp #sigma;template #sigma",
458 // 200, 0., .25, 200, 0., .25);
459 // AddHistogram(ID(kHistCmpSigmaTOFpSemi), "#sigma comparison;exp #sigma;template #sigma",
460 // 200, 0., .25, 200, 0., .25);
461 // AddHistogram(ID(kHistCmpSigmaTOFdSemi), "#sigma comparison;exp #sigma;template #sigma",
462 // 200, 0., .25, 200, 0., .25);
463
464 // Nsigma distributions
465 AddHistogram(ID(kHistNsigmaTPCTOF), "N_{#sigma,p} TPC-TOF;p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
466 100, 0., 10.,
467 100, -25., 25.,
468 200, -100., 100.);
469 AddHistogram(ID(kHistNsigmaTPCTOFPt), "N_{#sigma,p} TPC-TOF;p_{T} (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
470 100, 0., 10.,
471 100, -25., 25.,
472 200, -100., 100.);
473 AddHistogram(ID(kHistNsigmaTPCTOFUsed), "N_{#sigma,p} TPC-TOF;p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
474 100, 0., 10.,
475 100, -25., 25.,
476 200, -100., 100.);
477 AddHistogram(ID(kHistNsigmaTPCTOFUsedCentral), "N_{#sigma,p} TPC-TOF (central);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
478 100, 0., 10.,
479 100, -25., 25.,
480 200, -100., 100.);
481 AddHistogram(ID(kHistNsigmaTPCTOFUsedSemiCentral), "N_{#sigma,p} TPC-TOF (semi-central);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
482 100, 0., 10.,
483 100, -25., 25.,
484 200, -100., 100.);
485 AddHistogram(ID(kHistNsigmaTPCTOFUsedPt), "N_{#sigma,p} TPC-TOF;p_{T} (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
486 50, 0., 10.,
487 100, -25., 25.,
488 200, -100., 100.);
489 AddHistogram(ID(kHistNsigmaTPCTOFUsedPtCentral), "N_{#sigma,p} TPC-TOF;p_{T} (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
490 50, 0., 10.,
491 100, -25., 25.,
492 200, -100., 100.);
493 AddHistogram(ID(kHistNsigmaTPCTOFUsedPtSemiCentral), "N_{#sigma,p} TPC-TOF;p_{T} (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
494 50, 0., 10.,
495 100, -25., 25.,
496 200, -100., 100.);
497
498 AddHistogram(ID(kHistNsigmaTPCTOFUsedCentralMCe), "N_{#sigma,p} TPC-TOF (central, MC e);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
499 100, 0., 10.,
500 100, -25., 25.,
501 200, -100., 100.);
502 AddHistogram(ID(kHistNsigmaTPCTOFUsedCentralMCmu), "N_{#sigma,p} TPC-TOF (central, MC #mu);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
503 100, 0., 10.,
504 100, -25., 25.,
505 200, -100., 100.);
506 AddHistogram(ID(kHistNsigmaTPCTOFUsedCentralMCpi), "N_{#sigma,p} TPC-TOF (central, MC pi);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
507 100, 0., 10.,
508 100, -25., 25.,
509 200, -100., 100.);
510 AddHistogram(ID(kHistNsigmaTPCTOFUsedCentralMCk), "N_{#sigma,p} TPC-TOF (central, MC K);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
511 100, 0., 10.,
512 100, -25., 25.,
513 200, -100., 100.);
514 AddHistogram(ID(kHistNsigmaTPCTOFUsedCentralMCp), "N_{#sigma,p} TPC-TOF (central, MC p);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
515 100, 0., 10.,
516 100, -25., 25.,
517 200, -100., 100.);
518 AddHistogram(ID(kHistNsigmaTPCTOFUsedCentralMCd), "N_{#sigma,p} TPC-TOF (central, MC d);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
519 100, 0., 10.,
520 100, -25., 25.,
521 200, -100., 100.);
522
523 AddHistogram(ID(kHistNsigmaTPCTOFUsedSemiCentralMCe), "N_{#sigma,p} TPC-TOF (semi-central, MC e);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
524 100, 0., 10.,
525 100, -25., 25.,
526 200, -100., 100.);
527 AddHistogram(ID(kHistNsigmaTPCTOFUsedSemiCentralMCmu), "N_{#sigma,p} TPC-TOF (semi-central, MC #mu);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
528 100, 0., 10.,
529 100, -25., 25.,
530 200, -100., 100.);
531 AddHistogram(ID(kHistNsigmaTPCTOFUsedSemiCentralMCpi), "N_{#sigma,p} TPC-TOF (semi-central, MC pi);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
532 100, 0., 10.,
533 100, -25., 25.,
534 200, -100., 100.);
535 AddHistogram(ID(kHistNsigmaTPCTOFUsedSemiCentralMCk), "N_{#sigma,p} TPC-TOF (semi-central, MC K);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
536 100, 0., 10.,
537 100, -25., 25.,
538 200, -100., 100.);
539 AddHistogram(ID(kHistNsigmaTPCTOFUsedSemiCentralMCp), "N_{#sigma,p} TPC-TOF (semi-central, MC p);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
540 100, 0., 10.,
541 100, -25., 25.,
542 200, -100., 100.);
543 AddHistogram(ID(kHistNsigmaTPCTOFUsedSemiCentralMCd), "N_{#sigma,p} TPC-TOF (semi-central, MC d);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
544 100, 0., 10.,
545 100, -25., 25.,
546 200, -100., 100.);
547
548 AddHistogram(ID(kHistEvPlane), "default event plane;#Psi;counts",
549 100, -0. * TMath::Pi(), 1. * TMath::Pi());
550 AddHistogram(ID(kHistEvPlaneUsed), "default event plane;#Psi;counts",
551 100, -0. * TMath::Pi(), 1. * TMath::Pi(),
552 kClLast, -.5, kClLast-.5);
553 AddHistogram(ID(kHistEvPlaneCheck), "backup event plane;#Psi;counts",
554 100, -1. * TMath::Pi(), 1. * TMath::Pi());
555 AddHistogram(ID(kHistEvPlaneCheckUsed), "backup event plane;#Psi;counts",
556 100, -0. * TMath::Pi(), 1. * TMath::Pi(),
557 kClLast, -.5, kClLast-.5);
558 AddHistogram(ID(kHistEvPlane3), "3rd order event plane;#Psi;counts",
559 100, -0. * TMath::Pi(), 2./3. * TMath::Pi());
560 AddHistogram(ID(kHistEvPlaneCorr), "default - backup event plane;#Psi_{def};#Psi_{bak};event class",
561 100, -0. * TMath::Pi(), 1. * TMath::Pi(),
562 100, -0. * TMath::Pi(), 1. * TMath::Pi(),
563 kClLast, -.5, kClLast-.5);
564 AddHistogram(ID(kHistEvPlaneCorrNoTrgJets), "event plane w/ and w/o trigger jets;#Psi_{full};#Psi_{notrgjets};event class",
565 100, -0. * TMath::Pi(), 1. * TMath::Pi(),
566 100, -0. * TMath::Pi(), 1. * TMath::Pi(),
567 kClLast, -.5, kClLast-.5);
568 AddHistogram(ID(kHistEvPlaneCorrNoTrgJetsTrgd), "event plane w/ and w/o trigger jets;#Psi_{full};#Psi_{notrgjets};event class",
569 100, -0. * TMath::Pi(), 1. * TMath::Pi(),
570 100, -0. * TMath::Pi(), 1. * TMath::Pi(),
571 kClLast, -.5, kClLast-.5);
572
573 AddHistogram(ID(kHistJetPtCentral), "jet spectrum - central;p_{T}^{jet,ch} (GeV/#it{c});counts",
574 40, 0., 200.);
575 AddHistogram(ID(kHistJetPtSemi), "jet spectrum - semi-peripheral;p_{T}^{jet,ch} (GeV/#it{c});counts",
576 40, 0., 200.);
577
578 AddHistogram(ID(kHistEtaPhiTrgHad), "trg had;#varphi;#eta",
579 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
580 100, -2., 2.);
581 AddHistogram(ID(kHistEtaPhiTrgJet), "trg jet;#varphi;#eta",
582 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
583 100, -2., 2.);
584 AddHistogram(ID(kHistEtaPhiAssHad), "ass had;#varphi;#eta",
585 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
586 100, -2., 2.);
587 AddHistogram(ID(kHistEtaPhiAssProt), "ass proton;#varphi;#eta",
588 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
589 100, -2., 2.);
590
591 AddHistogram(ID(kHistPhiTrgJetEvPlane), "trg jet;#varphi - #Psi_{ev};centrality",
592 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
593 100, 0., 100.);
594 AddHistogram(ID(kHistPhiTrgHadEvPlane), "trg had;#varphi - #Psi_{ev};centrality",
595 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
596 100, 0., 100.);
597 AddHistogram(ID(kHistPhiRndTrgJetEvPlane), "rnd trg jet;#varphi - #Psi_{ev};centrality",
598 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
599 100, 0., 100.);
600 AddHistogram(ID(kHistPhiRndTrgHadEvPlane), "rnd trg had;#varphi - #Psi_{ev};centrality",
601 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
602 100, 0., 100.);
603 AddHistogram(ID(kHistPhiAssHadEvPlane), "ass had;#varphi - #Psi_{ev};centrality",
604 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
605 100, 0., 100.);
606 AddHistogram(ID(kHistPhiAssProtEvPlane), "ass prot;#varphi - #Psi_{ev};centrality",
607 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
608 100, 0., 100.);
609 AddHistogram(ID(kHistPhiAssHadVsEvPlane), "ass had;#Psi_{ev};#varphi;centrality",
610 100, -0. * TMath::Pi(), 1. * TMath::Pi(),
611 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
612 100, 0., 100.);
613
614 AddHistogram(ID(kHistPhiTrgJetEvPlane3), "trg jet;#varphi - #Psi_{ev};centrality",
615 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
616 100, 0., 100.);
617 AddHistogram(ID(kHistPhiTrgHadEvPlane3), "trg had;#varphi - #Psi_{ev};centrality",
618 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
619 100, 0., 100.);
620 AddHistogram(ID(kHistPhiAssHadEvPlane3), "ass had;#varphi - #Psi_{ev};centrality",
621 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
622 100, 0., 100.);
623 AddHistogram(ID(kHistPhiAssProtEvPlane3), "ass prot;#varphi - #Psi_{ev};centrality",
624 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
625 100, 0., 100.);
626
627 for (Int_t iCorr = 0; iCorr < kCorrLast; ++iCorr) {
628 for (Int_t iCl = 0; iCl < kClLast; ++iCl) {
629 for (Int_t iEv = 0; iEv < kEvLast; ++iEv) {
630 // we don't need the mixed event histograms for the embedded excess particles
631 if ((iCorr > kCorrRndJetProt) && (iEv == kEvMix))
632 continue;
633
634 GetHistCorr((CorrType_t) iCorr, (Class_t) iCl, (Ev_t) iEv) =
635 new AliHistCorr(Form("corr_%s_%s_%s", fkCorrTypeName[iCorr], fkClassName[iCl], fkEvName[iEv]), fOutputList);
636 }
637 }
638 }
639
640 PostData(1, fOutputList);
641}
642
643Bool_t AliAnalysisTaskJetProtonCorr::Notify()
644{
645 // actions to be taken upon notification about input file change
646
647 return AliAnalysisTaskSE::Notify();
648}
649
650void AliAnalysisTaskJetProtonCorr::SetParamsTOF()
651{
652 fParamsTOF =
653 dynamic_cast<AliTOFPIDParams*> (fOADBContainerTOF->GetObject(fRunNumber, "TOFparams"));
654
655 if (!fParamsTOF)
656 AliError(Form("failed to load TOF parameters for run %i", fRunNumber));
657 else if (fDebug > 2) {
658 printf("loaded TOF parameters for run %i\n",
659 fRunNumber);
660 printf(" intrinsic resolution: %f ps\n",
661 fParamsTOF->GetTOFresolution());
662 printf(" tail fraction: %f\n",
663 fParamsTOF->GetTOFtail());
664 printf(" start time method: %i\n",
665 fParamsTOF->GetStartTimeMethod());
666 printf(" TOF signal parametrization: %f, %f, %f, %f\n",
667 fParamsTOF->GetSigParams(0), fParamsTOF->GetSigParams(1),
668 fParamsTOF->GetSigParams(2), fParamsTOF->GetSigParams(3));
669 printf(" matching loss MC: %f%%\n",
670 fParamsTOF->GetTOFmatchingLossMC());
671 printf(" additional mismatch for MC: %f%%\n",
672 fParamsTOF->GetTOFadditionalMismForMC());
673 printf(" time offset: %f\n",
674 fParamsTOF->GetTOFtimeOffset());
675 }
676}
677
678void AliAnalysisTaskJetProtonCorr::UserExec(Option_t * /* option */)
679{
680 // actual work
681
682 // setup pointers to input data (null if unavailable)
683 // mcEvent: MC input
684 // esdEvent: ESD input
685 // outEvent: AOD output
686 // aodEvent: AOD input if available, otherwise AOD output
687
688 fMCEvent = this->MCEvent();
689 fESDEvent = dynamic_cast<AliESDEvent*>(this->InputEvent()); // could also be AOD input
690 AliAODEvent* outEvent = this->AODEvent();
691 fAODEvent = dynamic_cast<AliAODEvent*> (this->InputEvent());
692 if (!fAODEvent)
693 fAODEvent = outEvent;
694
695 if ((fDebug > 0) && fESDEvent)
696 printf("event: %s-%06i\n", CurrentFileName(), fESDEvent->GetEventNumberInFile());
697
698 // record number of sampled events and detect trigger contributions
699 FillH1(kHistStat, kStatSeen);
700 if (!DetectTriggers()) {
701 AliError("Failed to detect the triggers");
702 return;
703 }
704
705 if (!IsTrigger(kTriggerInt))
706 return;
707
708 FillH1(kHistStat, kStatTrg);
709
710 // prepare the event
711 // (make sure it is cleaned up in the end)
712 if (PrepareEvent()) {
713 FillH1(kHistStat, kStatUsed);
714 FillH1(kHistCentrality, fCentrality);
715 FillH1(kHistCentralityCheck, fCentralityCheck);
716
717 FillH1(kHistVertexZ, fZvtx);
718
719 // event cuts
720 if (TMath::Abs(fZvtx) > 10.)
721 goto stop;
722 if (GetCentrality() > 90.)
723 goto stop;
724
725 FillH1(kHistStat, kStatEvCuts);
726
727 // event category
728 DetectClasses();
729 if (IsClass(kClCentral))
730 FillH1(kHistStat, kStatCentral);
731 if (IsClass(kClSemiCentral))
732 FillH1(kHistStat, kStatSemiCentral);
733
734 FillH1(kHistEvPlane, fEventPlaneAngle);
735 FillH1(kHistEvPlaneCheck, fEventPlaneAngleCheck);
736 FillH1(kHistEvPlane3, fEventPlaneAngle3);
737 for (Int_t iClass = 0; iClass < kClLast; ++iClass) {
738 if (IsClass((Class_t) iClass)) {
739 FillH2(kHistCentralityUsed, fCentrality, iClass);
740 FillH2(kHistCentralityCheckUsed, fCentralityCheck, iClass);
741 FillH2(kHistEvPlaneUsed, fEventPlaneAngle, iClass);
742 FillH2(kHistEvPlaneCheckUsed, fEventPlaneAngleCheck, iClass);
743 FillH3(kHistEvPlaneCorr, fEventPlaneAngle, fEventPlaneAngleCheck, iClass);
744 }
745 }
746
747 Bool_t corrEta = fPIDResponse->UseTPCEtaCorrection();
748 Bool_t corrMult = fPIDResponse->UseTPCMultiplicityCorrection();
749
750 // select trigger particles and potential associated particles/protons
751 TObjArray trgArray[kTrgLast];
752 TObjArray assArray[kAssLast];
753
754 // prepare TOF response
755 TF1 fTOFsignal("fTOFsignal", &AliAnalysisTaskJetProtonCorr::TOFsignal, -2440., 2440., 4);
756 fTOFsignal.SetParameter(0, 1.);
757 fTOFsignal.SetParameter(1, 0.);
758 if (fParamsTOF) {
759 Float_t res = fParamsTOF->GetTOFresolution();
760 Float_t tail = fParamsTOF->GetTOFtail() * res;
761 fTOFsignal.SetParameter(2, res);
762 fTOFsignal.SetParameter(3, tail);
763 }
764 else {
765 fTOFsignal.SetParameter(2, 85.);
766 fTOFsignal.SetParameter(3, 80.);
767 }
768
769 // associate candidates
770 const Int_t nPrimTracksAss = fPrimTrackArrayAss ? fPrimTrackArrayAss->GetEntries() : 0;
771 FillH2(kHistCentralityVsMult, fCentrality, nPrimTracksAss);
772 for (Int_t iTrack = 0; iTrack < nPrimTracksAss; ++iTrack) {
773 AliVTrack *trk = (AliVTrack*) fPrimTrackArrayAss->At(iTrack);
774 FillH2(kHistSignalTPC, trk->P(), trk->GetTPCsignal());
775 // ??? pt or p?
776 FillH2(kHistSignalTOF, trk->P(), trk->GetTOFsignal() * 1.e-3); // ps -> ns
777 FillH3(kHistNsigmaTPCTOF,
778 trk->P(),
779 fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
780 fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
781 FillH3(kHistNsigmaTPCTOFPt,
782 trk->Pt(),
783 fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
784 fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
785
786 FillH3(kHistPhiAssHadVsEvPlane, fEventPlaneAngle, trk->Phi(), fCentrality);
787 Float_t phiRel = trk->Phi() - fEventPlaneAngle;
788 if (gRandom->Rndm() > .5)
789 phiRel -= TMath::Pi();
790 if (phiRel < 0.)
791 phiRel += 2. * TMath::Pi();
792 if (phiRel < 0.)
793 AliError(Form("phiRel = %f less than zero, from phi = %f, psi = %f",
794 phiRel, trk->Phi(), fEventPlaneAngle));
795 else if (phiRel > 2*TMath::Pi())
796 AliError(Form("phiRel = %f greater than 2pi, from phi = %f, psi = %f",
797 phiRel, trk->Phi(), fEventPlaneAngle));
798 Float_t phiRel3 = trk->Phi() - fEventPlaneAngle3;
799 if (phiRel3 < 0.)
800 phiRel3 += 2. * TMath::Pi();
801
802 if (AcceptAssoc(trk)) {
803 assArray[kAssHad].Add(trk);
804 FillH1(kHistEtaPhiAssHad, trk->Phi(), trk->Eta());
805 FillH2(kHistPhiAssHadEvPlane, phiRel, fCentrality);
806 FillH2(kHistPhiAssHadEvPlane3, phiRel3, fCentrality);
807 FillH3(kHistNsigmaTPCTOFUsed,
808 trk->P(),
809 fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
810 fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
811 FillH3(kHistNsigmaTPCTOFUsedPt,
812 trk->Pt(),
813 fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
814 fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
815 if (IsClass(kClCentral)) {
816 FillH3(kHistNsigmaTPCTOFUsedCentral,
817 trk->P(),
818 fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
819 fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
820 FillH3(kHistNsigmaTPCTOFUsedPtCentral,
821 trk->Pt(),
822 fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
823 fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
824 if (fMCEvent) {
825 if (AliMCParticle *mcPart = (AliMCParticle*) fMCEvent->GetTrack(TMath::Abs(trk->GetLabel()))) {
826 for (Int_t iParticle = 0; iParticle <= AliPID::kDeuteron; ++iParticle) {
827 if (TMath::Abs(mcPart->Particle()->GetPdgCode()) == AliPID::ParticleCode(iParticle))
828 FillH3(kHistNsigmaTPCTOFUsedCentralMCe,
829 trk->P(),
830 fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
831 fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton),
832 1., iParticle);
833 }
834 }
835 }
836 }
837 if (IsClass(kClSemiCentral)) {
838 FillH3(kHistNsigmaTPCTOFUsedSemiCentral,
839 trk->P(),
840 fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
841 fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
842 FillH3(kHistNsigmaTPCTOFUsedPtSemiCentral,
843 trk->Pt(),
844 fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
845 fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
846 if (fMCEvent) {
847 if (AliMCParticle *mcPart = (AliMCParticle*) fMCEvent->GetTrack(TMath::Abs(trk->GetLabel()))) {
848 for (Int_t iParticle = 0; iParticle <= AliPID::kDeuteron; ++iParticle) {
849 if (TMath::Abs(mcPart->Particle()->GetPdgCode()) == AliPID::ParticleCode(iParticle))
850 FillH3(kHistNsigmaTPCTOFUsedSemiCentralMCe,
851 trk->P(),
852 fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
853 fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton),
854 1., iParticle);
855 }
856 }
857 }
858 }
859 if (IsProton(trk)) {
860 assArray[kAssProt].Add(trk);
861 FillH1(kHistEtaPhiAssProt, trk->Phi(), trk->Eta());
862 FillH2(kHistPhiAssProtEvPlane, phiRel, fCentrality);
863 FillH2(kHistPhiAssProtEvPlane3, phiRel3, fCentrality);
864 }
865
866 // template generation
867 if (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC, trk) == AliPIDResponse::kDetPidOk) {
868 Double_t deltaTPC;
869 fPIDResponse->GetSignalDelta(AliPIDResponse::kTPC, trk, AliPID::kProton, deltaTPC);
870 if (IsClass(kClCentral))
871 FillH2(kHistDeltaTPC, trk->P(), deltaTPC);
872 if (IsClass(kClSemiCentral))
873 FillH2(kHistDeltaTPCSemi, trk->P(), deltaTPC);
874
875 Double_t expTPC = fPIDResponse->GetTPCResponse().GetExpectedSignal(trk, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, corrEta, corrMult);
876 Double_t expSigmaTPC = fPIDResponse->GetTPCResponse().GetExpectedSigma(trk, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, corrEta, corrMult);
877
878 // loop over particles
879 for (Int_t iParticle = 0; iParticle <= AliPID::kDeuteron; ++iParticle) {
880 Double_t expTPCx = fPIDResponse->GetTPCResponse().GetExpectedSignal(trk, AliPID::EParticleType(iParticle), AliTPCPIDResponse::kdEdxDefault, corrEta, corrMult);
881 Double_t expSigmaTPCx = fPIDResponse->GetTPCResponse().GetExpectedSigma(trk, AliPID::EParticleType(iParticle), AliTPCPIDResponse::kdEdxDefault, corrEta, corrMult);
882 Double_t rndTPCx = gRandom->Gaus(expTPCx, expSigmaTPCx);
883
884 if(IsClass(kClCentral))
885 FillH2(kHistNsigmaTPCe, trk->P(), (rndTPCx-expTPC)/expSigmaTPC, 1., iParticle);
886 if(IsClass(kClSemiCentral))
887 FillH2(kHistNsigmaTPCeSemi, trk->P(), (rndTPCx-expTPC)/expSigmaTPC, 1., iParticle);
888 }
889 }
890
891 if (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, trk) == AliPIDResponse::kDetPidOk) {
892 Double_t deltaTOF;
893 fPIDResponse->GetSignalDelta(AliPIDResponse::kTOF, trk, AliPID::kProton, deltaTOF);
894
895 AliTOFPIDResponse &tofResponse = fPIDResponse->GetTOFResponse();
896
897 Float_t p = trk->P();
898 if (const AliVVtrack *param = trk->GetInnerParam())
899 p = param->GetP();
900
901 Double_t expTOF = fPIDResponse->GetTOFResponse().GetExpectedSignal(trk, AliPID::kProton);
902 Double_t expSigmaTOF = fPIDResponse->GetTOFResponse().GetExpectedSigma(p, expTOF, AliPID::kProton);
903 Double_t length = trk->GetIntegratedLength() * 1.e-2; // cm -> m
904 Double_t tof = trk->GetTOFsignal() * 1.e-12; // ps -> s
905 Double_t beta = length / tof / TMath::C();
906
907 FillH2(kHistBetaTOF, p, beta);
908
909 // intrinsic TOF smearing
910 // Double_t signalSigma = TMath::Sqrt(fTOFsignal.Variance(-2000., 2000.));
911 Double_t signalSmear = fTOFsignal.GetRandom();
912 // t0 smearing
913 Float_t timezeroSigma = tofResponse.GetT0binRes(tofResponse.GetMomBin(p));
914 Double_t timezeroSmear = gRandom->Gaus(0., timezeroSigma);
915 // tracking smearing (default parameters, should be overwritten from OADB)
916 Double_t fPar[] = { 0.008, 0.008, 0.002, 40.0 };
917 if (fParamsTOF)
918 for (Int_t i = 0; i < 4; ++i)
919 fPar[i] = fParamsTOF->GetSigParams(i);
920
921 // loop over particles
922 for (Int_t iParticle = 0; iParticle <= AliPID::kDeuteron; ++iParticle) {
923 Double_t expTOFx = fPIDResponse->GetTOFResponse().GetExpectedSignal(trk, AliPID::EParticleType(iParticle));
924 // Double_t cmpSigmaTOFx = fPIDResponse->GetTOFResponse().GetExpectedSigma(p, expTOFx, AliPID::EParticleType(iParticle));
925
926 // tracking smearing
927 Double_t massx = AliPID::ParticleMassZ(AliPID::EParticleType(iParticle));
928 Double_t dppx = fPar[0] + fPar[1] * p + fPar[2] * massx / p;
929 Double_t expSigmaTOFx = dppx * expTOFx / (1.+ p * p / (massx * massx));
930 Double_t texpSmearx = gRandom->Gaus(0., TMath::Sqrt(expSigmaTOFx * expSigmaTOFx + fPar[3]*fPar[3]/p/p));
931 // Double_t tmpSigmaTOFx = TMath::Sqrt(signalSigma*signalSigma +
932 // timezeroSigma*timezeroSigma +
933 // expSigmaTOFx*expSigmaTOFx + fPar[3]*fPar[3]/p/p);
934 // printf("sigma comparison %i, %f: %f, %f, %f -> %f vs %f\n",
935 // iParticle, expTOFx,
936 // signalSigma, // signal
937 // timezeroSigma, // timezero
938 // expSigmaTOFx, // tracking
939 // tmpSigmaTOFx, // total
940 // cmpSigmaTOFx); // from PID response
941
942 // TOF signal
943 Double_t rndTOFx = expTOFx + signalSmear + timezeroSmear + texpSmearx;
944
945 if (IsClass(kClCentral)) {
946 FillH2(kHistNsigmaTOFe, trk->P(), (rndTOFx-expTOF)/expSigmaTOF, 1., iParticle);
947 FillH2(kHistDeltaTOFe, trk->P(), (rndTOFx-expTOF) * 1.e-3, 1., iParticle);
948 // FillH2(kHistExpSigmaTOFe, p, cmpSigmaTOFx * 1.e-3, 1., iParticle);
949 // FillH2(kHistCmpSigmaTOFe, cmpSigmaTOFx * 1.e-3, tmpSigmaTOFx * 1.e-3, 1., iParticle);
950 }
951 if (IsClass(kClSemiCentral)) {
952 FillH2(kHistNsigmaTOFeSemi, trk->P(), (rndTOFx-expTOF)/expSigmaTOF, 1., iParticle);
953 FillH2(kHistDeltaTOFeSemi, trk->P(), (rndTOFx-expTOF) * 1.e-3, 1., iParticle);
954 // FillH2(kHistCmpSigmaTOFeSemi, cmpSigmaTOFx * 1.e-3, tmpSigmaTOFx * 1.e-3, 1., iParticle);
955 // FillH2(kHistExpSigmaTOFeSemi, p, cmpSigmaTOFx * 1.e-3, 1., iParticle);
956 }
957 }
958
959 Double_t rndTOFmismatch = AliTOFPIDResponse::GetMismatchRandomValue(trk->Eta());
960
961 if(IsClass(kClCentral)) {
962 FillH2(kHistNsigmaTOFmismatch, trk->P(), (rndTOFmismatch - expTOF) / expSigmaTOF);
963 FillH2(kHistDeltaTOF, trk->P(), deltaTOF * 1.e-3); // ps -> ns
964 }
965 if(IsClass(kClSemiCentral)) {
966 FillH2(kHistNsigmaTOFmismatchSemi, trk->P(), (rndTOFmismatch - expTOF) / expSigmaTOF);
967 FillH2(kHistDeltaTOFSemi, trk->P(), deltaTOF * 1.e-3); // ps -> ns
968 }
969 }
970 }
971 }
972
973 const Int_t nPrimTracksTrg = fPrimTrackArrayTrg ? fPrimTrackArrayTrg->GetEntries() : 0;
974 for (Int_t iTrack = 0; iTrack < nPrimTracksTrg; ++iTrack) {
975 AliVTrack *trk = (AliVTrack*) fPrimTrackArrayTrg->At(iTrack);
976 if (AcceptTrigger(trk)) {
977 trgArray[kTrgHad].Add(trk);
978 FillH1(kHistEtaPhiTrgHad, trk->Phi(), trk->Eta());
979 }
980 }
981
982 // select trigger jet
983 // and remove them from the Q vector
984 Int_t nJets = fJetArray ? fJetArray->GetEntries() : 0;
985 const TVector2 *qVectorOrig = fEventplane->GetQVector();
986 TVector2 qVector;
987 if (qVectorOrig)
988 qVector = *qVectorOrig;
989
990 for (Int_t iJet = 0; iJet < nJets; ++iJet) {
991 AliAODJet *jet = (AliAODJet*) fJetArray->At(iJet);
992
993 if (AcceptTrigger(jet)) {
994 trgArray[kTrgJet].Add(jet);
995 FillH1(kHistEtaPhiTrgJet, jet->Phi(), jet->Eta());
996
997 if (qVectorOrig) {
998 Int_t nRefTracks = jet->GetRefTracks()->GetEntriesFast();
999 for (Int_t iTrack = 0; iTrack < nRefTracks; ++iTrack) {
1000 AliVTrack *track = (AliVTrack*) jet->GetRefTracks()->At(iTrack);
1001
1002 if (fEventplane) {
1003 TVector2 evplaneContrib(fEventplane->GetQContributionX(track),
1004 fEventplane->GetQContributionY(track));
1005 qVector -= evplaneContrib;
1006 }
1007 }
1008 }
1009 }
1010 }
1011 if (qVectorOrig) {
1012 for (Int_t iClass = 0; iClass < kClLast; ++iClass) {
1013 if (IsClass((Class_t) iClass)) {
1014 FillH3(kHistEvPlaneCorrNoTrgJets, qVectorOrig->Phi()/2., qVector.Phi()/2., iClass);
1015 if (trgArray[kTrgJet].GetEntriesFast() > 0)
1016 FillH3(kHistEvPlaneCorrNoTrgJetsTrgd, qVectorOrig->Phi()/2., qVector.Phi()/2., iClass);
1017 }
1018 }
1019 }
1020
1021 // invent a trigger jet/hadron and correlated associates
1022 Float_t pFraction = assArray[kAssHad].GetEntries() > 0 ?
1023 (Float_t) (assArray[kAssProt].GetEntries()) / assArray[kAssHad].GetEntries() :
1024 .5;
1025 GenerateRandom(&trgArray[kTrgJetRnd], &trgArray[kTrgHadRnd],
1026 &assArray[kAssHadJetExc], &assArray[kAssProtJetExc],
1027 &assArray[kAssHadHadExc], &assArray[kAssProtHadExc],
1028 pFraction);
1029
1030 // correlate, both same and mixed event
1031 for (Int_t iClass = 0; iClass < kClLast; ++iClass) {
1032 if (IsClass((Class_t) iClass)) {
1033 for (Int_t iTrg = 0; iTrg < kTrgLast; ++iTrg) {
1034 for (Int_t iAss = 0; iAss <= kAssProt; ++iAss) {
1035 // same event
1036 Correlate((Trg_t) iTrg, (Ass_t) iAss, (Class_t) iClass, kEvSame, &trgArray[iTrg], &assArray[iAss]);
1037
1038 // mixed event
1039 AliEventPool *pool = GetPool((Ass_t) iAss);
1040 if (pool && pool->IsReady()) {
1041 AliDebug(1, Form("----- using pool: %i %i %i -----", iClass, iTrg, iAss));
1042 const Int_t nEvents = pool->GetCurrentNEvents();
1043 for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
1044 TObjArray *assTracks = pool->GetEvent(iEvent);
1045 Correlate((Trg_t) iTrg, (Ass_t) iAss, (Class_t) iClass, kEvMix, &trgArray[iTrg], assTracks, 1./nEvents);
1046 }
1047 }
1048 }
1049 }
1050
1051 // fill event pool for mixing
1052 // >= 0: don't require a trigger in the event
1053 // >= 1: require a trigger in the event
1054 // if (trgArray[iTrg].GetEntries() >= 0) {
1055 for (Int_t iAss = 0; iAss <= kAssProt; ++iAss) {
1056 AliEventPool *pool = GetPool((Ass_t) iAss);
1057 if (pool) {
1058 pool->UpdatePool(CloneTracks(&assArray[iAss]));
1059 AliDebug(1, Form("----- updating pool: %i -----", iAss));
1060 if (fDebug > 0)
1061 pool->PrintInfo();
1062 }
1063 }
1064 // }
1065
1066 // correlate artificial triggers and associates
1067 Correlate(kCorrRndJetExcHad, (Class_t) iClass, kEvSame, &trgArray[kTrgJetRnd], &assArray[kAssHadJetExc]);
1068 Correlate(kCorrRndJetExcProt, (Class_t) iClass, kEvSame, &trgArray[kTrgJetRnd], &assArray[kAssProtJetExc]);
1069 Correlate(kCorrRndHadExcHad, (Class_t) iClass, kEvSame, &trgArray[kTrgHadRnd], &assArray[kAssHadHadExc]);
1070 Correlate(kCorrRndHadExcProt, (Class_t) iClass, kEvSame, &trgArray[kTrgHadRnd], &assArray[kAssProtHadExc]);
1071 }
1072 }
1073
1074 trgArray[kTrgJetRnd].Delete();
1075 trgArray[kTrgHadRnd].Clear();
1076 assArray[kAssHadJetExc].Delete();
1077 assArray[kAssProtJetExc].Clear();
1078 assArray[kAssHadHadExc].Delete();
1079 assArray[kAssProtHadExc].Clear();
1080 }
1081
1082 stop:
1083 CleanUpEvent();
1084
1085 PostData(1, fOutputList);
1086}
1087
1088void AliAnalysisTaskJetProtonCorr::Terminate(const Option_t * /* option */)
1089{
1090 // actions at task termination
1091
1092}
1093
1094void AliAnalysisTaskJetProtonCorr::PrintTask(Option_t *option, Int_t indent) const
1095{
1096 AliAnalysisTaskSE::PrintTask(option, indent);
1097
1098 std::cout << std::setw(indent) << " " << "using jet branch: " << fJetBranchName << std::endl;
1099}
1100
1101Bool_t AliAnalysisTaskJetProtonCorr::DetectTriggers()
1102{
1103 fTriggerMask = 0;
1104
1105 AliVEvent::EOfflineTriggerTypes physSel =
1106 (AliVEvent::EOfflineTriggerTypes) ((AliInputEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1107 // TString trgClasses = InputEvent()->GetFiredTriggerClasses();
1108
1109 // physics selection
1110 if (physSel & (AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral))
1111 MarkTrigger(kTriggerInt);
1112
1113 return kTRUE;
1114}
1115
1116Bool_t AliAnalysisTaskJetProtonCorr::DetectClasses()
1117{
1118 fClassMask = 0;
1119
1120 if (IsCentral())
1121 MarkClass(kClCentral);
1122
1123 if (IsSemiCentral())
1124 MarkClass(kClSemiCentral);
1125
1126 return kTRUE;
1127}
1128
1129Bool_t AliAnalysisTaskJetProtonCorr::PrepareEvent()
1130{
1131 Bool_t eventGood = kTRUE;
1132
1133 // check for run change
1134 if (fRunNumber != InputEvent()->GetRunNumber()) {
1135 fRunNumber = InputEvent()->GetRunNumber();
1136 SetParamsTOF();
1137 }
1138
1139 // retrieve z-vertex position
1140 fZvtx = 100.;
1141 const AliVVertex *vtx = InputEvent()->GetPrimaryVertex();
1142 if (vtx) {
1143 FillH1(kHistStat, kStatVtx);
1144 FillH1(kHistVertexNctb, vtx->GetNContributors());
1145 if (vtx->GetNContributors() >= 3.)
1146 fZvtx = vtx->GetZ();
1147 }
1148
1149 // retrieve centrality
1150 AliCentrality *eventCentrality = InputEvent()->GetCentrality();
1151 if (eventCentrality) {
1152 fCentrality = eventCentrality->GetCentralityPercentile("V0M");
1153 fCentralityCheck = eventCentrality->GetCentralityPercentile("TRK");
1154 if (fCentrality >= 0.) {
1155 FillH1(kHistStat, kStatCent);
1156 } else {
1157 // centrality estimation not reliable
1158 eventGood = kFALSE;
1159 fCentrality = 105.;
1160 }
1161 }
1162 else
1163 eventGood = kFALSE;
1164
1165 // retrieve event plane
1166 fEventplane = InputEvent()->GetEventplane();
1167 if (fEventplane) {
1168 fEventPlaneAngle3 = fEventplane->GetEventplane("V0", InputEvent(), 3);
1169
1170 if (fUseEvplaneV0) {
1171 fEventPlaneAngle = fEventplane->GetEventplane("V0", InputEvent());
1172 fEventPlaneAngleCheck = fEventplane->GetEventplane("Q");
1173 }
1174 else {
1175 fEventPlaneAngle = fEventplane->GetEventplane("Q");
1176 fEventPlaneAngleCheck = fEventplane->GetEventplane("V0", InputEvent());
1177 // use V0 event plane angle from flow task:
1178 // fEventPlaneAngleCheck = AliAnalysisTaskVnV0::GetPsi2V0A();
1179 // printf("V0A evplane = %f\n", fEventPlaneAngleCheck);
1180 // fEventPlaneAngleCheck = AliAnalysisTaskVnV0::GetPsi2V0C();
1181 // printf("V0C evplane = %f\n", fEventPlaneAngleCheck);
1182 }
1183
1184 // ensure angles to be in [0, ...)
1185 if (fEventPlaneAngle3 < 0)
1186 fEventPlaneAngle3 += 2.*TMath::Pi()/3.;
1187 if (fEventPlaneAngle < 0)
1188 fEventPlaneAngle += TMath::Pi();
1189 if (fEventPlaneAngleCheck < 0)
1190 fEventPlaneAngleCheck += TMath::Pi();
1191
1192 FillH1(kHistStat, kStatEvPlane);
1193 }
1194 else
1195 eventGood = kFALSE;
1196
1197 // retrieve PID
1198 fPIDResponse = ((AliInputEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
1199 if (fPIDResponse)
1200 FillH1(kHistStat, kStatPID);
1201 else
1202 eventGood = kFALSE;
1203
1204 // retrieve primary tracks
1205 if (fESDEvent) {
1206 fPrimTrackArrayAss = fCutsPrimAss->GetAcceptedTracks(fESDEvent);
1207 fPrimTrackArrayTrg = fCutsPrimTrg->GetAcceptedTracks(fESDEvent);
1208 if (fCutsPrimTrgConstrain) {
1209 TIter trkIter(fCutsPrimTrgConstrain->GetAcceptedTracks(fESDEvent));
1210 while (AliESDtrack *trk = (AliESDtrack*) trkIter()) {
1211 if (!fCutsPrimTrg->IsSelected(trk)) {
1212 AliESDtrack *track = (AliESDtrack*) fPrimConstrainedTrackArray->ConstructedAt(fPrimConstrainedTrackArray->GetEntriesFast());
1213 if(trk->GetConstrainedParam()) {
1214 track->Set(trk->GetConstrainedParam()->GetX(),
1215 trk->GetConstrainedParam()->GetAlpha(),
1216 trk->GetConstrainedParam()->GetParameter(),
1217 trk->GetConstrainedParam()->GetCovariance());
1218 }
1219 fPrimTrackArrayTrg->Add(track);
1220 }
1221 }
1222 }
1223 }
1224 else if (fAODEvent) {
1225 // associate candidates
1226 fPrimTrackArrayAss = new TObjArray();
1227 const Int_t nTracksAODAss = fAODEvent->GetNumberOfTracks();
1228 for (Int_t iTrack = 0; iTrack < nTracksAODAss; ++iTrack) {
1229 AliAODTrack *trk = fAODEvent->GetTrack(iTrack);
1230 if (trk->TestFilterBit(fAssFilterMask))
1231 fPrimTrackArrayAss->Add(trk);
1232 }
1233
1234 // trigger candidates
1235 fPrimTrackArrayTrg = new TObjArray();
1236 const Int_t nTracksAODTrg = fAODEvent->GetNumberOfTracks();
1237 for (Int_t iTrack = 0; iTrack < nTracksAODTrg; ++iTrack) {
1238 AliAODTrack *trk = fAODEvent->GetTrack(iTrack);
1239 if (trk->IsHybridGlobalConstrainedGlobal())
1240 fPrimTrackArrayTrg->Add(trk);
1241 }
1242 }
1243 else
1244 eventGood = kFALSE;
1245
1246 // retrieve jet array
1247 if (fAODEvent) {
1248 fJetArray = dynamic_cast<TClonesArray*> (fAODEvent->FindListObject(fJetBranchName));
1249 if (!fJetArray) {
1250 // still try the output event
1251 if (AODEvent())
1252 fJetArray = dynamic_cast<TClonesArray*> (AODEvent()->FindListObject(fJetBranchName));
1253 if (!fJetArray) {
1254 if (fErrorMsg > 0) {
1255 AliError(Form("no jet branch \"%s\" found, in the AODs are:", fJetBranchName));
1256 --fErrorMsg;
1257 }
1258 if (fDebug > 0) {
1259 fAODEvent->GetList()->Print();
1260 if (AODEvent())
1261 AODEvent()->GetList()->Print();
1262 }
1263 }
1264 }
1265 }
1266
1267 // retrieve event pool for the event category
1268 if (eventGood) {
1269 for (Int_t iAss = 0; iAss <= kAssProt; ++iAss) {
1270 AliEventPoolManager *mgr = GetPoolMgr((Ass_t) iAss);
1271 GetPool((Ass_t) iAss) =
1272 // mgr ? mgr->GetEventPool(fCentrality, fZvtx, fEventPlaneAngle) : 0x0;
1273 mgr ? mgr->GetEventPool(fCentrality, fZvtx) : 0x0;
1274 }
1275 }
1276
1277 return eventGood;
1278}
1279
1280Bool_t AliAnalysisTaskJetProtonCorr::AcceptTrigger(const AliVTrack *trg)
1281{
1282 if (TMath::Abs(trg->Eta()) > fHadEtaMax)
1283 return kFALSE;
1284
1285 // check pt interval
1286 Bool_t acceptPt =
1287 (trg->Pt() >= fTrgPartPtMin) &&
1288 (trg->Pt() <= fTrgPartPtMax);
1289
1290 // for semi-central events check phi w.r.t. event plane
1291 Bool_t acceptOrientation = IsSemiCentral() ?
1292 AcceptAngleToEvPlane(trg->Phi(), GetEventPlaneAngle()) :
1293 kTRUE;
1294
1295 if (acceptPt) {
1296 Float_t phiRel = trg->Phi() - fEventPlaneAngle;
1297 if (gRandom->Rndm() > .5)
1298 phiRel -= TMath::Pi();
1299 if (phiRel < 0.)
1300 phiRel += 2. * TMath::Pi();
1301 Float_t phiRel3 = trg->Phi() - fEventPlaneAngle3;
1302 if (phiRel3 < 0.)
1303 phiRel3 += 2. * TMath::Pi();
1304 FillH2(kHistPhiTrgHadEvPlane, phiRel, fCentrality);
1305 FillH2(kHistPhiTrgHadEvPlane3, phiRel3, fCentrality);
1306 }
1307
1308 return (acceptPt && acceptOrientation);
1309}
1310
1311AliVTrack* AliAnalysisTaskJetProtonCorr::GetLeadingTrack(const AliAODJet *jet) const
1312{
1313 // return leading track within a jet
1314
1315 // check contributing tracks
1316 Int_t nJetTracks = jet->GetRefTracks()->GetEntriesFast();
1317 Int_t iLeadingTrack = -1;
1318 Float_t ptLeadingTrack = 0.;
1319 for (Int_t iTrack=0; iTrack < nJetTracks; ++iTrack) {
1320 AliAODTrack *track = (AliAODTrack*) jet->GetRefTracks()->At(iTrack);
1321 // find the leading track
1322 if (track->Pt() > ptLeadingTrack) {
1323 ptLeadingTrack = track->Pt();
1324 iLeadingTrack = iTrack;
1325 }
1326 }
1327
1328 // retrieve the leading track
1329 return (AliAODTrack*) jet->GetRefTracks()->At(iLeadingTrack);
1330}
1331
1332Bool_t AliAnalysisTaskJetProtonCorr::AcceptTrigger(const AliAODJet *trg)
1333{
1334 // restrict eta
1335 if (TMath::Abs(trg->Eta()) > fTrgJetEtaMax)
1336 return kFALSE;
1337
1338 // reject too small jets
1339 if (trg->EffectiveAreaCharged() < fTrgJetAreaMin)
1340 return kFALSE;
1341
1342 // require hard leading track
1343 // (leading track biased jet sample)
1344 // but reject jets with too high pt constituents
1345 const Float_t ptLeadTrack = GetLeadingTrack(trg)->Pt();
1346 if ((ptLeadTrack < fTrgJetLeadTrkPtMin) ||
1347 (ptLeadTrack > fTrgJetLeadTrkPtMax))
1348 return kFALSE;
1349
1350 // check for jet orientation w.r.t. event plane
1351 // (constrained only for semi-central events)
1352 Bool_t acceptOrientation = IsSemiCentral() ?
1353 AcceptAngleToEvPlane(trg->Phi(), GetEventPlaneAngle()) :
1354 kTRUE;
1355
1356 // check for jet pt
1357 Bool_t acceptPt =
1358 (trg->Pt() >= fTrgJetPtMin) &&
1359 (trg->Pt() <= fTrgJetPtMax);
1360
1361 if (acceptPt) {
1362 // store the azimuthal distribution relative to the event plane
1363 // for the jets in the pt range of interest
1364 Float_t phiRel = trg->Phi() - fEventPlaneAngle;
1365 if (gRandom->Rndm() > .5)
1366 phiRel -= TMath::Pi();
1367 if (phiRel < 0.)
1368 phiRel += 2. * TMath::Pi();
1369 Float_t phiRel3 = trg->Phi() - fEventPlaneAngle3;
1370 if (phiRel3 < 0.)
1371 phiRel3 += 2. * TMath::Pi();
1372 FillH2(kHistPhiTrgJetEvPlane, phiRel, fCentrality);
1373 FillH2(kHistPhiTrgJetEvPlane3, phiRel3, fCentrality);
1374 }
1375
1376 if (acceptOrientation) {
1377 // spectrum for cross checks
1378 if (IsClass(kClCentral))
1379 FillH1(kHistJetPtCentral, trg->Pt());
1380 if (IsClass(kClSemiCentral))
1381 FillH1(kHistJetPtSemi, trg->Pt());
1382 }
1383
1384 return (acceptPt && acceptOrientation);
1385}
1386
1387Bool_t AliAnalysisTaskJetProtonCorr::AcceptAssoc(const AliVTrack *trk) const
1388{
1389 if ((trk->Pt() > fAssPartPtMin) && (trk->Pt() < fAssPartPtMax) && (TMath::Abs(trk->Eta()) < fHadEtaMax))
1390 if (fRequirePID) {
1391 if ((fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC, trk) == AliPIDResponse::kDetPidOk) &&
1392 (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, trk) == AliPIDResponse::kDetPidOk) &&
1393 (trk->GetTPCsignalN() >= 60.) && (trk->GetTPCsignal() >= 10))
1394 return kTRUE;
1395 else
1396 return kFALSE;
1397 }
1398 else
1399 return kTRUE;
1400 else
1401 return kFALSE;
1402}
1403
1404Bool_t AliAnalysisTaskJetProtonCorr::IsProton(const AliVTrack *trk) const
1405{
1406 if ((fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC, trk) != AliPIDResponse::kDetPidOk) ||
1407 (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, trk) != AliPIDResponse::kDetPidOk) ||
1408 (trk->GetTPCsignalN() < 60.) || (trk->GetTPCsignal() < 10))
1409 return kFALSE;
1410
1411 Double_t nSigmaProtonTPC = fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton);
1412 Double_t nSigmaProtonTOF = fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton);
1413
1414 if ((TMath::Abs(nSigmaProtonTPC) <= 2.) && (TMath::Abs(nSigmaProtonTOF) <= 2.)) {
1415 return kTRUE;
1416 }
1417
1418 return kFALSE;
1419}
1420
1421Bool_t AliAnalysisTaskJetProtonCorr::AcceptAngleToEvPlane(Float_t phi, Float_t psi) const
1422{
1423 Float_t deltaPhi = phi - psi;
1424 while (deltaPhi < 0.)
1425 deltaPhi += 2 * TMath::Pi();
1426
1427 // map to interval [-pi/2, pi/2)
1428 deltaPhi = std::fmod(deltaPhi + TMath::Pi()/2., TMath::Pi()) - TMath::Pi()/2.;
1429 // printf("delta phi = %f with phi = %f, psi = %f\n", deltaPhi, phi, psi);
1430
1431 if (TMath::Abs(deltaPhi) < fTrgAngleToEvPlane)
1432 return kTRUE;
1433 else
1434 return kFALSE;
1435}
1436
1437Float_t AliAnalysisTaskJetProtonCorr::GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1,
1438 Float_t phi2, Float_t pt2, Float_t charge2,
1439 Float_t radius, Float_t bSign) const
1440{
1441 // calculates dphistar
1442 // from AliUEHistograms
1443
1444 Float_t dphistar = phi1 - phi2
1445 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1)
1446 + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
1447
1448 static const Double_t kPi = TMath::Pi();
1449
1450 // circularity
1451 if (dphistar > kPi)
1452 dphistar = kPi * 2 - dphistar;
1453 if (dphistar < -kPi)
1454 dphistar = -kPi * 2 - dphistar;
1455 if (dphistar > kPi) // might look funny but is needed
1456 dphistar = kPi * 2 - dphistar;
1457
1458 return dphistar;
1459}
1460
1461
1462Bool_t AliAnalysisTaskJetProtonCorr::AcceptTwoTracks(const AliVParticle *trgPart, const AliVParticle *assPart) const
1463{
1464 // apply two track pair cut
1465
1466 Float_t phi1 = trgPart->Phi();
1467 Float_t pt1 = trgPart->Pt();
1468 Float_t charge1 = trgPart->Charge();
1469
1470 Float_t phi2 = assPart->Phi();
1471 Float_t pt2 = assPart->Pt();
1472 Float_t charge2 = assPart->Charge();
1473
1474 Float_t deta = trgPart->Eta() - assPart->Eta();
1475
1476 Float_t bSign = (InputEvent()->GetMagneticField() > 0) ? 1 : -1;
1477
1478 // optimization
1479 if (TMath::Abs(deta) < fCutsTwoTrackEff) {
1480 // check first boundaries to see if is worth to loop and find the minimum
1481 Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 0.8, bSign);
1482 Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 2.5, bSign);
1483
1484 const Float_t kLimit = fCutsTwoTrackEff * 3;
1485
1486 Float_t dphistarminabs = 1e5;
1487 // Float_t dphistarmin = 1e5;
1488 if ((TMath::Abs(dphistar1) < kLimit) ||
1489 (TMath::Abs(dphistar2) < kLimit) ||
1490 (dphistar1 * dphistar2 < 0)) {
1491 for (Double_t rad=0.8; rad<2.51; rad+=0.01) {
1492 Float_t dphistarabs = TMath::Abs(GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign));
1493 if (dphistarabs < dphistarminabs)
1494 dphistarminabs = dphistarabs;
1495 }
1496
1497 if (dphistarminabs < fCutsTwoTrackEff)
1498 return kFALSE;
1499 }
1500 }
1501
1502 return kTRUE;
1503}
1504
1505Bool_t AliAnalysisTaskJetProtonCorr::Correlate(CorrType_t corr, Class_t cl, Ev_t ev,
1506 TCollection *trgArray, TCollection *assArray, Float_t weight)
1507{
1508 AliHistCorr *histCorr = GetHistCorr(corr, cl, ev);
1509 if (!histCorr) {
1510 AliError(Form("no correlation histograms for corr %i, cl %i, ev %i",
1511 corr, cl, ev));
1512 return kFALSE;
1513 }
1514
1515 TIter assIter(assArray);
1516 while (TObject *assObj = assIter()) {
1517 if (dynamic_cast<AliVParticle*> (assObj)) {
1518 AliVParticle *assPart = (AliVParticle*) assObj;
1519 histCorr->Ass(assPart->Phi(), assPart->Eta(), weight);
1520 }
1521 else if (dynamic_cast<TLorentzVector*> (assObj)) {
1522 TLorentzVector *assVec = (TLorentzVector*) assObj;
1523 histCorr->Ass(assVec->Phi(), assVec->Eta(), weight);
1524 }
1525 }
1526
1527 TIter trgIter(trgArray);
1528 while (TObject *trgObj = trgIter()) {
1529 if (AliVParticle *trgPart = dynamic_cast<AliVParticle*> (trgObj)) {
1530 // count the trigger
1531 histCorr->Trigger(trgPart->Phi(), trgPart->Eta(), weight);
1532
1533 // loop over associates
1534 assIter.Reset();
1535 while (TObject *assObj = assIter()) {
1536 if (AliVParticle *assPart = dynamic_cast<AliVParticle*> (assObj)) {
1537 if (AcceptTwoTracks(trgPart, assPart))
1538 histCorr->Fill(trgPart, assPart, weight);
1539 }
1540 else if (TLorentzVector *assVec = dynamic_cast<TLorentzVector*> (assObj)) {
1541 AliFatal(Form("got %p, but not implemented", assVec));
1542 }
1543 }
1544 }
1545 else if (TLorentzVector *trgVec = dynamic_cast<TLorentzVector*> (trgObj)) {
1546 // count the trigger
1547 histCorr->Trigger(trgVec->Phi(), trgVec->Eta(), weight);
1548
1549 // loop over associates
1550 assIter.Reset();
1551 while (TObject *assObj = assIter()) {
1552 if (AliVParticle *assPart = dynamic_cast<AliVParticle*> (assObj)) {
1553 histCorr->Fill(trgVec, assPart, weight);
1554 }
1555 else if (TLorentzVector *assVec = dynamic_cast<TLorentzVector*> (assObj)) {
1556 histCorr->Fill(trgVec, assVec, weight);
1557 }
1558 }
1559 }
1560 }
1561
1562 return kTRUE;
1563}
1564
1565Bool_t AliAnalysisTaskJetProtonCorr::Correlate(Trg_t trg, Ass_t ass, Class_t cl, Ev_t ev,
1566 TCollection *trgArray, TCollection *assArray, Float_t weight)
1567{
1568 if ((trg < 0) || (trg > kTrgJetRnd))
1569 AliFatal(Form("wrong request for correlation with trigger: %d", trg));
1570
1571 if ((ass < 0) || (ass > kAssProt))
1572 AliFatal(Form("wrong request for correlation with associate: %d", ass));
1573
1574 CorrType_t corr = (CorrType_t) (2 * trg + ass);
1575
1576 return Correlate(corr, cl, ev, trgArray, assArray, weight);
1577}
1578
1579Bool_t AliAnalysisTaskJetProtonCorr::GenerateRandom(TCollection *trgJetArray,
1580 TCollection *trgHadArray,
1581 TCollection *assHadJetArray,
1582 TCollection *assProtJetArray,
1583 TCollection *assHadHadArray,
1584 TCollection *assProtHadArray,
1585 Float_t pFraction)
1586{
1587 // generate random direction
1588
1589 const Int_t nPart = gRandom->Poisson(fToyMeanNoPart);
1590
1591 Float_t trgJetEta = gRandom->Uniform(-fTrgJetEtaMax, fTrgJetEtaMax);
1592 Float_t trgJetPhi = 0.;
1593 do {
1594 trgJetPhi = IsClass(kClSemiCentral) ?
1595 fTrgJetPhiModSemi->GetRandom() :
1596 fTrgJetPhiModCent->GetRandom();
1597
1598 trgJetPhi += fEventPlaneAngle;
1599 if (trgJetPhi < 0.)
1600 trgJetPhi += 2. * TMath::Pi();
1601 else if (trgJetPhi > 2. * TMath::Pi())
1602 trgJetPhi -= 2. * TMath::Pi();
1603
1604 Float_t phiRel = trgJetPhi- fEventPlaneAngle;
1605 if (phiRel < 0.)
1606 phiRel += 2. * TMath::Pi();
1607 FillH2(kHistPhiRndTrgJetEvPlane, phiRel, fCentrality);
1608 } while (IsClass(kClSemiCentral) &&
1609 !AcceptAngleToEvPlane(trgJetPhi, GetEventPlaneAngle()));
1610
1611 // generate one trigger particle
1612 TLorentzVector *trgJet = new TLorentzVector();
1613 trgJet->SetPtEtaPhiM(50., trgJetEta, trgJetPhi, .14);
1614 trgJetArray->Add(trgJet);
1615
1616 // generate direction for away side
1617 Float_t trgJetEtaAway = gRandom->Uniform(-fHadEtaMax, fHadEtaMax);
1618 Float_t trgJetPhiAway = trgJetPhi + TMath::Pi() + gRandom->Gaus(0., fToySmearPhi);
1619
1620 // generate associated particles
1621 // with proton/hadron ratio observed in this event
1622 for (Int_t iPart = 0; iPart < nPart; ++iPart) {
1623 Float_t deltaEta, deltaPhi;
1624 const Float_t r = fToyRadius;
1625 do {
1626 deltaEta = gRandom->Uniform(-r, r);
1627 deltaPhi = gRandom->Uniform(-r, r);
1628 } while ((deltaEta * deltaEta + deltaPhi * deltaPhi) > (r * r));
1629
1630 TLorentzVector *assPart = new TLorentzVector();
1631 Float_t eta = trgJetEtaAway + deltaEta;
1632 Float_t phi = trgJetPhiAway + deltaPhi;
1633 if (TMath::Abs(eta) > fHadEtaMax) {
1634 delete assPart;
1635 continue;
1636 }
1637 while (phi < 0.)
1638 phi += 2. * TMath::Pi();
1639 while (phi > 2. * TMath::Pi())
1640 phi -= 2. * TMath::Pi();
1641 assPart->SetPtEtaPhiM(3., eta, phi, .14);
1642 assHadJetArray->Add(assPart);
1643 if (gRandom->Uniform() < pFraction)
1644 assProtJetArray->Add(assPart);
1645 }
1646
1647 // trigger hadron
1648 Float_t trgHadEta = gRandom->Uniform(-fHadEtaMax, fHadEtaMax);
1649 Float_t trgHadPhi = 0.;
1650 do {
1651 trgHadPhi = IsClass(kClSemiCentral) ?
1652 fTrgHadPhiModSemi->GetRandom() :
1653 fTrgHadPhiModCent->GetRandom();
1654
1655 trgHadPhi += fEventPlaneAngle;
1656 if (trgHadPhi < 0.)
1657 trgHadPhi += 2. * TMath::Pi();
1658 else if (trgHadPhi > 2. * TMath::Pi())
1659 trgHadPhi -= 2. * TMath::Pi();
1660
1661 Float_t phiRel = trgHadPhi - fEventPlaneAngle;
1662 if (phiRel < 0.)
1663 phiRel += 2. * TMath::Pi();
1664 FillH2(kHistPhiRndTrgHadEvPlane, phiRel, fCentrality);
1665 } while (IsClass(kClSemiCentral) &&
1666 !AcceptAngleToEvPlane(trgHadPhi, GetEventPlaneAngle()));
1667
1668 // generate one trigger particle
1669 TLorentzVector *trgHad = new TLorentzVector();
1670 trgHad->SetPtEtaPhiM(50., trgHadEta, trgHadPhi, .14);
1671 trgHadArray->Add(trgHad);
1672
1673 // generate direction for away side
1674 Float_t trgHadEtaAway = gRandom->Uniform(-fHadEtaMax, fHadEtaMax);
1675 Float_t trgHadPhiAway = trgHadPhi + TMath::Pi() + gRandom->Gaus(0., fToySmearPhi);
1676
1677 // generate associated particles
1678 // with proton/hadron ratio observed in this event
1679 for (Int_t iPart = 0; iPart < nPart; ++iPart) {
1680 Float_t deltaEta, deltaPhi;
1681 const Float_t r = fToyRadius;
1682 do {
1683 deltaEta = gRandom->Uniform(-r, r);
1684 deltaPhi = gRandom->Uniform(-r, r);
1685 } while ((deltaEta * deltaEta + deltaPhi * deltaPhi) > (r * r));
1686
1687 TLorentzVector *assPart = new TLorentzVector();
1688 Float_t eta = trgHadEtaAway + deltaEta;
1689 Float_t phi = trgHadPhiAway + deltaPhi;
1690 if (TMath::Abs(eta) > fHadEtaMax) {
1691 delete assPart;
1692 continue;
1693 }
1694 while (phi < 0.)
1695 phi += 2. * TMath::Pi();
1696 while (phi > 2. * TMath::Pi())
1697 phi -= 2. * TMath::Pi();
1698 assPart->SetPtEtaPhiM(3., eta, phi, .14);
1699 assHadHadArray->Add(assPart);
1700 if (gRandom->Uniform() < pFraction)
1701 assProtHadArray->Add(assPart);
1702 }
1703
1704 return kTRUE;
1705}
1706
1707Bool_t AliAnalysisTaskJetProtonCorr::CleanUpEvent()
1708{
1709 if (fAODEvent) {
1710 delete fPrimTrackArrayAss;
1711 fPrimTrackArrayAss = 0x0;
1712 delete fPrimTrackArrayTrg;
1713 fPrimTrackArrayTrg = 0x0;
1714 }
1715
1716 fPrimConstrainedTrackArray->Delete();
1717
1718 return kTRUE;
1719}
1720
1721TObjArray* AliAnalysisTaskJetProtonCorr::CloneTracks(TObjArray* tracks) const
1722{
1723 TObjArray* tracksClone = new TObjArray;
1724 tracksClone->SetOwner(kTRUE);
1725
1726 Int_t nTracks = tracks->GetEntriesFast();
1727 for (Int_t i = 0; i < nTracks; i++) {
1728 // WARNING: TObject::Clone() is very!!! expensive, unusable
1729 // tracksClone->Add(particle->Clone());
1730
1731 // if (AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (tracks->At(i)))
1732 // tracksClone->Add(new AliESDtrack(*esdTrack));
1733 // else if (AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*> (tracks->At(i)))
1734 // tracksClone->Add(new AliAODTrack(*aodTrack));
1735
1736 if (const AliVTrack *track = dynamic_cast<const AliVTrack*> (tracks->At(i)))
1737 tracksClone->Add(new AliPartCorr(*track));
1738 }
1739
1740 return tracksClone;
1741}
1742
1743AliAnalysisTaskJetProtonCorr::AliHistCorr::AliHistCorr(TString name, TList *outputList) :
1744 TNamed(name, name),
1745 fOutputList(outputList),
1746 fHistStat(0x0),
1747 fHistCorrPhi(0x0),
1748 fHistCorrPhi2(0x0),
1749 fHistCorrEtaPhi(0x0),
1750 fHistCorrAvgEtaPhi(0x0),
1751 fHistCorrTrgEtaPhi(0x0),
1752 fHistCorrAssEtaPhi(0x0),
1753 fHistDphiLo(-TMath::Pi() / 2.),
1754 fHistDphiNbins(120),
1755 fHistDetaNbins(120)
1756{
1757 // ctor
1758
1759 fHistStat = new TH1F(Form("%s_stat", name.Data()), "statistics",
1760 1, .5, 1.5);
1761
1762 fHistCorrPhi = new TH1F(Form("%s_phi", name.Data()), ";#Delta#phi",
1763 fHistDphiNbins, fHistDphiLo, 2.*TMath::Pi() + fHistDphiLo);
1764 fHistCorrPhi->Sumw2();
1765 // fHistCorrPhi2 = new TH2F(Form("%s_phi2", name.Data()), ";#phi_{trg};#phi_{ass}",
1766 // 120, 0.*TMath::Pi(), 2.*TMath::Pi(),
1767 // 120, 0.*TMath::Pi(), 2.*TMath::Pi());
1768 // fHistCorrPhi2->Sumw2();
1769 fHistCorrEtaPhi = new TH2F(Form("%s_etaphi", name.Data()), ";#Delta#phi;#Delta#eta",
1770 fHistDphiNbins, fHistDphiLo, 2*TMath::Pi() + fHistDphiLo,
1771 fHistDetaNbins, -2., 2.);
1772 fHistCorrEtaPhi->Sumw2();
1773 // fHistCorrAvgEtaPhi = new TH2F(Form("%s_avgetaphi", name.Data()), ";#Delta#phi;avg #eta",
1774 // fHistDphiNbins, fHistDphiLo, 2*TMath::Pi() + fHistDphiLo,
1775 // fHistDetaNbins, -2., 2.);
1776 // fHistCorrAvgEtaPhi->Sumw2();
1777 // fHistCorrTrgEtaPhi = new TH2F(Form("%s_trg_etaphi", name.Data()), ";#varphi;#eta",
1778 // 120, 0., 2*TMath::Pi(),
1779 // 120, -1., 1.);
1780 // fHistCorrTrgEtaPhi->Sumw2();
1781 // fHistCorrAssEtaPhi = new TH2F(Form("%s_ass_etaphi", name.Data()), ";#varphi;#eta",
1782 // 120, 0., 2*TMath::Pi(),
1783 // 120, -1., 1.);
1784 // fHistCorrAssEtaPhi->Sumw2();
1785
1786 fOutputList->Add(fHistStat);
1787 fOutputList->Add(fHistCorrPhi);
1788 // fOutputList->Add(fHistCorrPhi2);
1789 fOutputList->Add(fHistCorrEtaPhi);
1790 // fOutputList->Add(fHistCorrAvgEtaPhi);
1791 // fOutputList->Add(fHistCorrTrgEtaPhi);
1792 // fOutputList->Add(fHistCorrAssEtaPhi);
1793}
1794
1795AliAnalysisTaskJetProtonCorr::AliHistCorr::~AliHistCorr()
1796{
1797 // dtor
1798}
1799
1800void AliAnalysisTaskJetProtonCorr::AliHistCorr::Fill(AliVParticle *trgPart, AliVParticle *assPart, Float_t weight)
1801{
1802 // fill correlation histograms from trigger particle and associate particle
1803
1804 Float_t deltaEta = assPart->Eta() - trgPart->Eta();
1805 Float_t avgEta = (assPart->Eta() + trgPart->Eta()) / 2.;
1806 Float_t deltaPhi = assPart->Phi() - trgPart->Phi();
1807 while (deltaPhi > (2.*TMath::Pi() + fHistDphiLo))
1808 deltaPhi -= 2. * TMath::Pi();
1809 while (deltaPhi < fHistDphiLo)
1810 deltaPhi += 2. * TMath::Pi();
1811
1812 fHistCorrPhi->Fill(deltaPhi, weight);
1813 if (fHistCorrPhi2)
1814 fHistCorrPhi2->Fill(trgPart->Phi(), assPart->Phi(), weight);
1815 fHistCorrEtaPhi->Fill(deltaPhi, deltaEta, weight);
1816 if (fHistCorrAvgEtaPhi)
1817 fHistCorrAvgEtaPhi->Fill(deltaPhi, avgEta, weight);
1818}
1819
1820void AliAnalysisTaskJetProtonCorr::AliHistCorr::Fill(TLorentzVector *trgPart, AliVParticle *assPart, Float_t weight)
1821{
1822 // fill correlation histograms from trigger direction and associate particle
1823
1824 Float_t deltaEta = assPart->Eta() - trgPart->Eta();
1825 Float_t avgEta = (assPart->Eta() + trgPart->Eta()) / 2.;
1826 Float_t deltaPhi = assPart->Phi() - trgPart->Phi();
1827 while (deltaPhi > (2.*TMath::Pi() + fHistDphiLo))
1828 deltaPhi -= 2. * TMath::Pi();
1829 while (deltaPhi < fHistDphiLo)
1830 deltaPhi += 2. * TMath::Pi();
1831
1832 fHistCorrPhi->Fill(deltaPhi, weight);
1833 if (fHistCorrPhi2) {
1834 Float_t trgPhi = trgPart->Phi();
1835 if (trgPhi < 0)
1836 trgPhi += 2. * TMath::Pi();
1837 fHistCorrPhi2->Fill(trgPhi, assPart->Phi(), weight);
1838 }
1839 fHistCorrEtaPhi->Fill(deltaPhi, deltaEta, weight);
1840 if (fHistCorrAvgEtaPhi)
1841 fHistCorrAvgEtaPhi->Fill(deltaPhi, avgEta, weight);
1842}
1843
1844void AliAnalysisTaskJetProtonCorr::AliHistCorr::Fill(TLorentzVector *trgPart, TLorentzVector *assPart, Float_t weight)
1845{
1846 // fill correlation histograms from trigger direction and associate direction
1847
1848 Float_t deltaEta = assPart->Eta() - trgPart->Eta();
1849 Float_t avgEta = (assPart->Eta() + trgPart->Eta()) / 2.;
1850 Float_t deltaPhi = assPart->Phi() - trgPart->Phi();
1851 if (deltaPhi > (2.*TMath::Pi() + fHistDphiLo))
1852 deltaPhi -= 2. * TMath::Pi();
1853 else if (deltaPhi < fHistDphiLo)
1854 deltaPhi += 2. * TMath::Pi();
1855
1856 fHistCorrPhi->Fill(deltaPhi, weight);
1857 if (fHistCorrPhi2) {
1858 Float_t trgPhi = trgPart->Phi();
1859 if (trgPhi < 0)
1860 trgPhi += 2. * TMath::Pi();
1861 Float_t assPhi = assPart->Phi();
1862 if (assPhi < 0)
1863 assPhi += 2. * TMath::Pi();
1864 fHistCorrPhi2->Fill(trgPhi, assPhi, weight);
1865 }
1866 fHistCorrEtaPhi->Fill(deltaPhi, deltaEta, weight);
1867 if (fHistCorrAvgEtaPhi)
1868 fHistCorrAvgEtaPhi->Fill(deltaPhi, avgEta, weight);
1869}
1870
1871// ----- histogram management -----
1872TH1* AliAnalysisTaskJetProtonCorr::AddHistogram(Hist_t hist, const char *hid, TString title,
1873 Int_t xbins, Float_t xmin, Float_t xmax,
1874 Int_t binType)
1875{
1876 TString hName;
1877 hName.Form("%s_%s", fShortTaskId, hid);
1878 hName.ToLower();
1879 TH1 *h = 0x0;
1880 if (binType == 0)
1881 h = new TH1I(hName.Data(), title,
1882 xbins, xmin, xmax);
1883 else
1884 h = new TH1F(hName.Data(), title,
1885 xbins, xmin, xmax);
1886 GetHistogram(hist) = h;
1887 fOutputList->Add(h);
1888 return h;
1889}
1890
1891TH2* AliAnalysisTaskJetProtonCorr::AddHistogram(Hist_t hist, const char *hid, TString title,
1892 Int_t xbins, Float_t xmin, Float_t xmax,
1893 Int_t ybins, Float_t ymin, Float_t ymax,
1894 Int_t binType)
1895{
1896 TString hName;
1897 hName.Form("%s_%s", fShortTaskId, hid);
1898 hName.ToLower();
1899 TH1 *h = 0x0;
1900 if (binType == 0)
1901 h = GetHistogram(hist) = new TH2I(hName.Data(), title,
1902 xbins, xmin, xmax,
1903 ybins, ymin, ymax);
1904 else
1905 h = GetHistogram(hist) = new TH2F(hName.Data(), title,
1906 xbins, xmin, xmax,
1907 ybins, ymin, ymax);
1908 fOutputList->Add(h);
1909 return (TH2*) h;
1910}
1911
1912TH3* AliAnalysisTaskJetProtonCorr::AddHistogram(Hist_t hist, const char *hid, TString title,
1913 Int_t xbins, Float_t xmin, Float_t xmax,
1914 Int_t ybins, Float_t ymin, Float_t ymax,
1915 Int_t zbins, Float_t zmin, Float_t zmax,
1916 Int_t binType)
1917{
1918 TString hName;
1919 hName.Form("%s_%s", fShortTaskId, hid);
1920 hName.ToLower();
1921 TH1 *h = 0x0;
1922 if (binType == 0)
1923 h = GetHistogram(hist) = new TH3I(hName.Data(), title,
1924 xbins, xmin, xmax,
1925 ybins, ymin, ymax,
1926 zbins, zmin, zmax);
1927 else
1928 h = GetHistogram(hist) = new TH3F(hName.Data(), title,
1929 xbins, xmin, xmax,
1930 ybins, ymin, ymax,
1931 zbins, zmin, zmax);
1932 fOutputList->Add(h);
1933 return (TH3*) h;
1934}