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