add common abstract interface classes for flat and fat ESDs
[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"
23536a91 10
11// analysis framework
12#include "AliAnalysisManager.h"
13#include "AliInputEventHandler.h"
14#include "AliVEvent.h"
ebdf1dbb 15#include "AliVTrack.h"
23536a91 16#include "AliVTrdTrack.h"
17#include "AliVVertex.h"
18#include "AliPIDResponse.h"
19#include "AliEventPoolManager.h"
ebdf1dbb 20#include "AliOADBContainer.h"
21#include "AliTOFPIDParams.h"
22#include "AliAnalysisTaskVnV0.h"
23536a91 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
3082e8ab 48#include <iostream>
49#include <cmath>
50
23536a91 51AliAnalysisTaskJetProtonCorr::AliAnalysisTaskJetProtonCorr(const char *name) :
52 AliAnalysisTaskSE(name),
53 fMCEvent(0x0),
54 fESDEvent(0x0),
55 fAODEvent(0x0),
ebdf1dbb 56 fRunNumber(-1),
57 fOADBContainerTOF(0x0),
58 fParamsTOF(0x0),
59 fEventplane(0x0),
23536a91 60 fTriggerMask(0),
61 fClassMask(0),
62 fCentrality(100.),
63 fCentralityCheck(100.),
64 fZvtx(0.),
65 fPIDResponse(0x0),
ebdf1dbb 66 fEventPlaneAngle(5.),
67 fEventPlaneAngleCheck(5.),
a4c7ef28 68 fEventPlaneAngle3(5.),
69 fPrimTrackArrayAss(0x0),
70 fPrimTrackArrayTrg(0x0),
71 fPrimConstrainedTrackArray(new TClonesArray("AliESDtrack", 100)),
23536a91 72 fJetArray(0x0),
73 fPoolMgr(),
74 fPool(),
75 fHistCorr(0x0),
d57e2346 76 fErrorMsg(10),
23536a91 77 fOutputList(),
78 fHist(),
79 fShortTaskId("jet_prot_corr"),
3082e8ab 80 fUseStandardCuts(kTRUE),
a4c7ef28 81 fUseEvplaneV0(kFALSE),
82 fCutsPrimTrg(0x0),
83 fCutsPrimTrgConstrain(new AliESDtrackCuts()),
84 fCutsPrimAss(0x0),
3082e8ab 85 fCutsTwoTrackEff(0.02),
d57e2346 86 fAssFilterMask(1 << 10),
687dfd6b 87 fRequirePID(kTRUE),
d57e2346 88 fTrgJetEtaMax(0.45),
89 fHadEtaMax(0.8),
23536a91 90 fTrgPartPtMin(6.),
91 fTrgPartPtMax(8.),
92 fTrgJetPtMin(50.),
93 fTrgJetPtMax(80.),
d57e2346 94 fTrgJetLeadTrkPtMin(6.),
95 fTrgJetLeadTrkPtMax(100.),
96 fTrgJetAreaMin(0.6 * TMath::Pi() * 0.2*0.2),
23536a91 97 fAssPartPtMin(2.),
3082e8ab 98 fAssPartPtMax(4.),
ebdf1dbb 99 fTrgAngleToEvPlane(TMath::Pi() / 4.),
687dfd6b 100 fToyMeanNoPart(1.),
101 fToyRadius(.8),
102 fToySmearPhi(.2),
ebdf1dbb 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()))
23536a91 107{
108 // default ctor
109
110 fkCorrTypeName[kCorrHadHad] = "hh";
111 fkCorrTypeName[kCorrHadProt] = "hp";
112 fkCorrTypeName[kCorrJetHad] = "jh";
113 fkCorrTypeName[kCorrJetProt] = "jp";
ebdf1dbb 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";
23536a91 122
123 fkClassName[kClCentral] = "cent";
124 fkClassName[kClSemiCentral] = "semi";
2751fb9a 125 // fkClassName[kClDijet] = "dijet";
23536a91 126
127 fkEvName[kEvSame] = "same";
128 fkEvName[kEvMix] = "mixed";
129
a4c7ef28 130 // track cuts for associates
3082e8ab 131 if (fUseStandardCuts) {
a4c7ef28 132 fCutsPrimAss = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
3082e8ab 133 } else {
a4c7ef28 134 fCutsPrimAss = new AliESDtrackCuts();
3082e8ab 135
136 // this is taken from PWGJE track cuts
137 TFormula *f1NClustersTPCLinearPtDep = new TFormula("f1NClustersTPCLinearPtDep","70.+30./20.*x");
a4c7ef28 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);
3082e8ab 145 // ITS
a4c7ef28 146 fCutsPrimAss->SetRequireITSRefit(kTRUE);
3082e8ab 147 //accept secondaries
a4c7ef28 148 fCutsPrimAss->SetMaxDCAToVertexXY(2.4);
149 fCutsPrimAss->SetMaxDCAToVertexZ(3.2);
150 fCutsPrimAss->SetDCAToVertex2D(kTRUE);
3082e8ab 151 //reject fakes
a4c7ef28 152 fCutsPrimAss->SetMaxChi2PerClusterITS(36);
153 fCutsPrimAss->SetMaxChi2TPCConstrainedGlobal(36);
3082e8ab 154
a4c7ef28 155 fCutsPrimAss->SetRequireSigmaToVertex(kFALSE);
156
157 fCutsPrimAss->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
3082e8ab 158 }
23536a91 159
d57e2346 160 fCutsPrimAss->SetEtaRange(-fHadEtaMax, fHadEtaMax);
a4c7ef28 161 fCutsPrimAss->SetPtRange(0.15, 1E+15);
162
163 // track cuts for triggers
164 fCutsPrimTrg = new AliESDtrackCuts(*fCutsPrimAss);
23536a91 165
a4c7ef28 166 // azimuthal modulation for triggers
d57e2346 167 fTrgJetPhiModCent->SetParameter(0, .10);
168 fTrgJetPhiModSemi->SetParameter(0, .10);
a4c7ef28 169 fTrgHadPhiModCent->SetParameter(0, .04);
170 fTrgHadPhiModSemi->SetParameter(0, .10);
ebdf1dbb 171
23536a91 172 // event mixing pool
3082e8ab 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 };
23536a91 178 Int_t nCentralityBins = sizeof(centralityBins)/sizeof(centralityBins[0]);
179
3082e8ab 180 Double_t vertexBins[] = {
181 -10., -8., -6., -4., -2., 0., 2., 4., 6., 8., 10.
182 };
23536a91 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
2751fb9a 190 const Int_t poolSize = 10; // unused by current event pool implementation
a4c7ef28 191 const Int_t trackDepth = 10000; // number of tracks to maintain in mixing buffer
2751fb9a 192 const Float_t trackDepthFraction = 0.1;
193 const Int_t targetEvents = 1;
579a71d2 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 }
23536a91 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
ebdf1dbb 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
23536a91 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));
ebdf1dbb 243 histStat->GetXaxis()->SetBinLabel(ID(kStatVtx));
23536a91 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));
3082e8ab 249 histStat->GetXaxis()->SetBinLabel(ID(kStatCentral));
250 histStat->GetXaxis()->SetBinLabel(ID(kStatSemiCentral));
23536a91 251
ebdf1dbb 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
23536a91 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));
2751fb9a 264 // hist->GetYaxis()->SetBinLabel(LAB(kClDijet));
23536a91 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));
2751fb9a 272 // hist->GetYaxis()->SetBinLabel(LAB(kClDijet));
ebdf1dbb 273 AddHistogram(ID(kHistCentralityVsMult), "centrality - multiplicity;centrality percentile (%);N_{prim}",
274 100, 0., 100.,
275 100, 0., 2500.);
23536a91 276
0ba09863 277 AddHistogram(ID(kHistSignalTPC), "TPC dE/dx;p (GeV/#it{c});dE/dx (arb. units)",
3082e8ab 278 100, 0., 10., 200, 0., 300.);
0ba09863 279 AddHistogram(ID(kHistSignalTOF), "TOF time of flight;p_{T} (GeV/#it{c});t (ns)",
3082e8ab 280 100, 0., 10., 200, 0., 50.);
0ba09863 281 AddHistogram(ID(kHistBetaTOF), "TOF beta;p (GeV/#it{c}); #beta",
3082e8ab 282 100, 0., 10.,
283 100, 0., 1.);
ebdf1dbb 284
0ba09863 285 AddHistogram(ID(kHistDeltaTPC), "TPC dE/dx;p (GeV/#it{c});dE/dx (arb. units)",
3082e8ab 286 100, 0., 10., 200, -100., 100.);
0ba09863 287 AddHistogram(ID(kHistDeltaTPCSemi), "TPC dE/dx;p (GeV/#it{c});dE/dx (arb. units)",
3082e8ab 288 100, 0., 10., 200, -100., 100.);
0ba09863 289 AddHistogram(ID(kHistDeltaTOF), "TOF time of flight;p (GeV/#it{c});t (ns)",
3082e8ab 290 100, 0., 10., 200, -2., 2.);
0ba09863 291 AddHistogram(ID(kHistDeltaTOFSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
3082e8ab 292 100, 0., 10., 200, -2., 2.);
293
ebdf1dbb 294 // Nsigma templates - central
0ba09863 295 AddHistogram(ID(kHistNsigmaTPCe), "TPC N_{#sigma,p} - e hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 296 100, 0., 10.,
297 100, -25., 25.);
0ba09863 298 AddHistogram(ID(kHistNsigmaTPCmu), "TPC N_{#sigma,p} - #mu hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 299 100, 0., 10.,
300 100, -25., 25.);
0ba09863 301 AddHistogram(ID(kHistNsigmaTPCpi), "TPC N_{#sigma,p} - #pi hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 302 100, 0., 10.,
303 100, -25., 25.);
0ba09863 304 AddHistogram(ID(kHistNsigmaTPCk), "TPC N_{#sigma,p} - K hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 305 100, 0., 10.,
306 100, -25., 25.);
0ba09863 307 AddHistogram(ID(kHistNsigmaTPCp), "TPC N_{#sigma,p} - p hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 308 100, 0., 10.,
309 100, -25., 25.);
0ba09863 310 AddHistogram(ID(kHistNsigmaTPCd), "TPC N_{#sigma,p} - d hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 311 100, 0., 10.,
312 100, -25., 25.);
0ba09863 313 AddHistogram(ID(kHistNsigmaTPCe_e), "TPC N_{#sigma,p} - e hypothesis (id. e);p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 314 100, 0., 10.,
315 100, -25., 25.);
316
0ba09863 317 AddHistogram(ID(kHistNsigmaTOFe), "TOF N_{#sigma,p} - e hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 318 100, 0., 10.,
ebdf1dbb 319 200, -100., 100.);
0ba09863 320 AddHistogram(ID(kHistNsigmaTOFmu), "TOF N_{#sigma,p} - #mu hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 321 100, 0., 10.,
ebdf1dbb 322 200, -100., 100.);
0ba09863 323 AddHistogram(ID(kHistNsigmaTOFpi), "TOF N_{#sigma,p} - #pi hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 324 100, 0., 10.,
ebdf1dbb 325 200, -100., 100.);
0ba09863 326 AddHistogram(ID(kHistNsigmaTOFk), "TOF N_{#sigma,p} - K hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 327 100, 0., 10.,
ebdf1dbb 328 200, -100., 100.);
0ba09863 329 AddHistogram(ID(kHistNsigmaTOFp), "TOF N_{#sigma,p} - p hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 330 100, 0., 10.,
ebdf1dbb 331 200, -100., 100.);
0ba09863 332 AddHistogram(ID(kHistNsigmaTOFd), "TOF N_{#sigma,p} - d hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 333 100, 0., 10.,
ebdf1dbb 334 200, -100., 100.);
0ba09863 335 AddHistogram(ID(kHistNsigmaTOFmismatch), "TOF N_{#sigma,p} - mismatch;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 336 100, 0., 10.,
ebdf1dbb 337 200, -100., 100.);
3082e8ab 338
ebdf1dbb 339 // Nsigma templates - semi-central
0ba09863 340 AddHistogram(ID(kHistNsigmaTPCeSemi), "TPC N_{#sigma,p} - e hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 341 100, 0., 10.,
342 100, -25., 25.);
0ba09863 343 AddHistogram(ID(kHistNsigmaTPCmuSemi), "TPC N_{#sigma,p} - #mu hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 344 100, 0., 10.,
345 100, -25., 25.);
0ba09863 346 AddHistogram(ID(kHistNsigmaTPCpiSemi), "TPC N_{#sigma,p} - #pi hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 347 100, 0., 10.,
348 100, -25., 25.);
0ba09863 349 AddHistogram(ID(kHistNsigmaTPCkSemi), "TPC N_{#sigma,p} - K hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 350 100, 0., 10.,
351 100, -25., 25.);
0ba09863 352 AddHistogram(ID(kHistNsigmaTPCpSemi), "TPC N_{#sigma,p} - p hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 353 100, 0., 10.,
354 100, -25., 25.);
0ba09863 355 AddHistogram(ID(kHistNsigmaTPCdSemi), "TPC N_{#sigma,p} - d hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 356 100, 0., 10.,
357 100, -25., 25.);
0ba09863 358 AddHistogram(ID(kHistNsigmaTPCe_eSemi), "TPC N_{#sigma,p} - e hypothesis (id. e);p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 359 100, 0., 10.,
360 100, -25., 25.);
361
0ba09863 362 AddHistogram(ID(kHistNsigmaTOFeSemi), "TOF N_{#sigma,p} - e hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 363 100, 0., 10.,
ebdf1dbb 364 200, -100., 100.);
0ba09863 365 AddHistogram(ID(kHistNsigmaTOFmuSemi), "TOF N_{#sigma,p} - #mu hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 366 100, 0., 10.,
ebdf1dbb 367 200, -100., 100.);
0ba09863 368 AddHistogram(ID(kHistNsigmaTOFpiSemi), "TOF N_{#sigma,p} - #pi hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 369 100, 0., 10.,
ebdf1dbb 370 200, -100., 100.);
0ba09863 371 AddHistogram(ID(kHistNsigmaTOFkSemi), "TOF N_{#sigma,p} - K hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 372 100, 0., 10.,
ebdf1dbb 373 200, -100., 100.);
0ba09863 374 AddHistogram(ID(kHistNsigmaTOFpSemi), "TOF N_{#sigma,p} - p hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 375 100, 0., 10.,
ebdf1dbb 376 200, -100., 100.);
0ba09863 377 AddHistogram(ID(kHistNsigmaTOFdSemi), "TOF N_{#sigma,p} - d hypothesis;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 378 100, 0., 10.,
ebdf1dbb 379 200, -100., 100.);
0ba09863 380 AddHistogram(ID(kHistNsigmaTOFmismatchSemi), "TOF N_{#sigma,p} - mismatch;p (GeV/#it{c});N_{#sigma,p}",
3082e8ab 381 100, 0., 10.,
ebdf1dbb 382 200, -100., 100.);
3082e8ab 383
384 // delta templates
0ba09863 385 AddHistogram(ID(kHistDeltaTOFe), "TOF #Delta;p (GeV/#it{c});t (ns)",
3082e8ab 386 100, 0., 10., 200, -2., 2.);
0ba09863 387 AddHistogram(ID(kHistDeltaTOFmu), "TOF #Delta;p (GeV/#it{c});t (ns)",
3082e8ab 388 100, 0., 10., 200, -2., 2.);
0ba09863 389 AddHistogram(ID(kHistDeltaTOFpi), "TOF #Delta;p (GeV/#it{c});t (ns)",
3082e8ab 390 100, 0., 10., 200, -2., 2.);
0ba09863 391 AddHistogram(ID(kHistDeltaTOFk), "TOF #Delta;p (GeV/#it{c});t (ns)",
3082e8ab 392 100, 0., 10., 200, -2., 2.);
0ba09863 393 AddHistogram(ID(kHistDeltaTOFp), "TOF #Delta;p (GeV/#it{c});t (ns)",
3082e8ab 394 100, 0., 10., 200, -2., 2.);
0ba09863 395 AddHistogram(ID(kHistDeltaTOFd), "TOF #Delta;p (GeV/#it{c});t (ns)",
3082e8ab 396 100, 0., 10., 200, -2., 2.);
397
0ba09863 398 AddHistogram(ID(kHistDeltaTOFeSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
3082e8ab 399 100, 0., 10., 200, -2., 2.);
0ba09863 400 AddHistogram(ID(kHistDeltaTOFmuSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
3082e8ab 401 100, 0., 10., 200, -2., 2.);
0ba09863 402 AddHistogram(ID(kHistDeltaTOFpiSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
3082e8ab 403 100, 0., 10., 200, -2., 2.);
0ba09863 404 AddHistogram(ID(kHistDeltaTOFkSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
3082e8ab 405 100, 0., 10., 200, -2., 2.);
0ba09863 406 AddHistogram(ID(kHistDeltaTOFpSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
3082e8ab 407 100, 0., 10., 200, -2., 2.);
0ba09863 408 AddHistogram(ID(kHistDeltaTOFdSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
3082e8ab 409 100, 0., 10., 200, -2., 2.);
410
411 // sigma comparisons
0ba09863 412 // AddHistogram(ID(kHistExpSigmaTOFe), "TOF time of flight;p (GeV/#it{c});t (ns)",
ebdf1dbb 413 // 100, 0., 10., 200, 0., .25);
0ba09863 414 // AddHistogram(ID(kHistExpSigmaTOFmu), "TOF time of flight;p (GeV/#it{c});t (ns)",
ebdf1dbb 415 // 100, 0., 10., 200, 0., .25);
0ba09863 416 // AddHistogram(ID(kHistExpSigmaTOFpi), "TOF time of flight;p (GeV/#it{c});t (ns)",
ebdf1dbb 417 // 100, 0., 10., 200, 0., .25);
0ba09863 418 // AddHistogram(ID(kHistExpSigmaTOFk), "TOF time of flight;p (GeV/#it{c});t (ns)",
ebdf1dbb 419 // 100, 0., 10., 200, 0., .25);
0ba09863 420 // AddHistogram(ID(kHistExpSigmaTOFp), "TOF time of flight;p (GeV/#it{c});t (ns)",
ebdf1dbb 421 // 100, 0., 10., 200, 0., .25);
0ba09863 422 // AddHistogram(ID(kHistExpSigmaTOFd), "TOF time of flight;p (GeV/#it{c});t (ns)",
ebdf1dbb 423 // 100, 0., 10., 200, 0., .25);
424
0ba09863 425 // AddHistogram(ID(kHistExpSigmaTOFeSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
ebdf1dbb 426 // 100, 0., 10., 200, 0., .25);
0ba09863 427 // AddHistogram(ID(kHistExpSigmaTOFmuSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
ebdf1dbb 428 // 100, 0., 10., 200, 0., .25);
0ba09863 429 // AddHistogram(ID(kHistExpSigmaTOFpiSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
ebdf1dbb 430 // 100, 0., 10., 200, 0., .25);
0ba09863 431 // AddHistogram(ID(kHistExpSigmaTOFkSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
ebdf1dbb 432 // 100, 0., 10., 200, 0., .25);
0ba09863 433 // AddHistogram(ID(kHistExpSigmaTOFpSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
ebdf1dbb 434 // 100, 0., 10., 200, 0., .25);
0ba09863 435 // AddHistogram(ID(kHistExpSigmaTOFdSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
ebdf1dbb 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);
3082e8ab 463
464 // Nsigma distributions
0ba09863 465 AddHistogram(ID(kHistNsigmaTPCTOF), "N_{#sigma,p} TPC-TOF;p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
3082e8ab 466 100, 0., 10.,
467 100, -25., 25.,
ebdf1dbb 468 200, -100., 100.);
0ba09863 469 AddHistogram(ID(kHistNsigmaTPCTOFPt), "N_{#sigma,p} TPC-TOF;p_{T} (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
3082e8ab 470 100, 0., 10.,
471 100, -25., 25.,
ebdf1dbb 472 200, -100., 100.);
0ba09863 473 AddHistogram(ID(kHistNsigmaTPCTOFUsed), "N_{#sigma,p} TPC-TOF;p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
3082e8ab 474 100, 0., 10.,
475 100, -25., 25.,
ebdf1dbb 476 200, -100., 100.);
0ba09863 477 AddHistogram(ID(kHistNsigmaTPCTOFUsedCentral), "N_{#sigma,p} TPC-TOF (central);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
3082e8ab 478 100, 0., 10.,
479 100, -25., 25.,
ebdf1dbb 480 200, -100., 100.);
0ba09863 481 AddHistogram(ID(kHistNsigmaTPCTOFUsedSemiCentral), "N_{#sigma,p} TPC-TOF (semi-central);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
3082e8ab 482 100, 0., 10.,
483 100, -25., 25.,
ebdf1dbb 484 200, -100., 100.);
0ba09863 485 AddHistogram(ID(kHistNsigmaTPCTOFUsedPt), "N_{#sigma,p} TPC-TOF;p_{T} (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
3082e8ab 486 50, 0., 10.,
487 100, -25., 25.,
ebdf1dbb 488 200, -100., 100.);
0ba09863 489 AddHistogram(ID(kHistNsigmaTPCTOFUsedPtCentral), "N_{#sigma,p} TPC-TOF;p_{T} (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
3082e8ab 490 50, 0., 10.,
491 100, -25., 25.,
ebdf1dbb 492 200, -100., 100.);
0ba09863 493 AddHistogram(ID(kHistNsigmaTPCTOFUsedPtSemiCentral), "N_{#sigma,p} TPC-TOF;p_{T} (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
3082e8ab 494 50, 0., 10.,
495 100, -25., 25.,
ebdf1dbb 496 200, -100., 100.);
3082e8ab 497
579a71d2 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
3082e8ab 548 AddHistogram(ID(kHistEvPlane), "default event plane;#Psi;counts",
23536a91 549 100, -0. * TMath::Pi(), 1. * TMath::Pi());
3082e8ab 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);
a4c7ef28 558 AddHistogram(ID(kHistEvPlane3), "3rd order event plane;#Psi;counts",
559 100, -0. * TMath::Pi(), 2./3. * TMath::Pi());
ebdf1dbb 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",
3082e8ab 569 100, -0. * TMath::Pi(), 1. * TMath::Pi(),
23536a91 570 100, -0. * TMath::Pi(), 1. * TMath::Pi(),
571 kClLast, -.5, kClLast-.5);
572
579a71d2 573 AddHistogram(ID(kHistJetPtCentral), "jet spectrum - central;p_{T}^{jet,ch} (GeV/#it{c});counts",
3082e8ab 574 40, 0., 200.);
579a71d2 575 AddHistogram(ID(kHistJetPtSemi), "jet spectrum - semi-peripheral;p_{T}^{jet,ch} (GeV/#it{c});counts",
23536a91 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
ebdf1dbb 591 AddHistogram(ID(kHistPhiTrgJetEvPlane), "trg jet;#varphi - #Psi_{ev};centrality",
592 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
21be932f 593 100, 0., 100.);
ebdf1dbb 594 AddHistogram(ID(kHistPhiTrgHadEvPlane), "trg had;#varphi - #Psi_{ev};centrality",
595 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
21be932f 596 100, 0., 100.);
687dfd6b 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.);
ebdf1dbb 603 AddHistogram(ID(kHistPhiAssHadEvPlane), "ass had;#varphi - #Psi_{ev};centrality",
604 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
21be932f 605 100, 0., 100.);
ebdf1dbb 606 AddHistogram(ID(kHistPhiAssProtEvPlane), "ass prot;#varphi - #Psi_{ev};centrality",
607 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
21be932f 608 100, 0., 100.);
ebdf1dbb 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(),
21be932f 612 100, 0., 100.);
ebdf1dbb 613
a4c7ef28 614 AddHistogram(ID(kHistPhiTrgJetEvPlane3), "trg jet;#varphi - #Psi_{ev};centrality",
615 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
21be932f 616 100, 0., 100.);
a4c7ef28 617 AddHistogram(ID(kHistPhiTrgHadEvPlane3), "trg had;#varphi - #Psi_{ev};centrality",
618 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
21be932f 619 100, 0., 100.);
a4c7ef28 620 AddHistogram(ID(kHistPhiAssHadEvPlane3), "ass had;#varphi - #Psi_{ev};centrality",
621 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
21be932f 622 100, 0., 100.);
a4c7ef28 623 AddHistogram(ID(kHistPhiAssProtEvPlane3), "ass prot;#varphi - #Psi_{ev};centrality",
624 100, -0. * TMath::Pi(), 2. * TMath::Pi(),
21be932f 625 100, 0., 100.);
a4c7ef28 626
23536a91 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) {
2751fb9a 630 // we don't need the mixed event histograms for the embedded excess particles
3fa79056 631 if ((iCorr > kCorrRndJetProt) && (iEv == kEvMix))
2751fb9a 632 continue;
633
23536a91 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
ebdf1dbb 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
23536a91 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
ebdf1dbb 717 FillH1(kHistVertexZ, fZvtx);
718
23536a91 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();
3082e8ab 729 if (IsClass(kClCentral))
730 FillH1(kHistStat, kStatCentral);
731 if (IsClass(kClSemiCentral))
732 FillH1(kHistStat, kStatSemiCentral);
23536a91 733
ebdf1dbb 734 FillH1(kHistEvPlane, fEventPlaneAngle);
735 FillH1(kHistEvPlaneCheck, fEventPlaneAngleCheck);
a4c7ef28 736 FillH1(kHistEvPlane3, fEventPlaneAngle3);
23536a91 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);
ebdf1dbb 741 FillH2(kHistEvPlaneUsed, fEventPlaneAngle, iClass);
742 FillH2(kHistEvPlaneCheckUsed, fEventPlaneAngleCheck, iClass);
743 FillH3(kHistEvPlaneCorr, fEventPlaneAngle, fEventPlaneAngleCheck, iClass);
23536a91 744 }
745 }
746
3082e8ab 747 Bool_t corrEta = fPIDResponse->UseTPCEtaCorrection();
748 Bool_t corrMult = fPIDResponse->UseTPCMultiplicityCorrection();
749
23536a91 750 // select trigger particles and potential associated particles/protons
751 TObjArray trgArray[kTrgLast];
752 TObjArray assArray[kAssLast];
753
579a71d2 754 // prepare TOF response
3082e8ab 755 TF1 fTOFsignal("fTOFsignal", &AliAnalysisTaskJetProtonCorr::TOFsignal, -2440., 2440., 4);
756 fTOFsignal.SetParameter(0, 1.);
757 fTOFsignal.SetParameter(1, 0.);
ebdf1dbb 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 }
3082e8ab 768
a4c7ef28 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);
ebdf1dbb 774 FillH2(kHistSignalTPC, trk->P(), trk->GetTPCsignal());
775 // ??? pt or p?
776 FillH2(kHistSignalTOF, trk->P(), trk->GetTOFsignal() * 1.e-3); // ps -> ns
23536a91 777 FillH3(kHistNsigmaTPCTOF,
3082e8ab 778 trk->P(),
779 fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
780 fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
781 FillH3(kHistNsigmaTPCTOFPt,
23536a91 782 trk->Pt(),
783 fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
784 fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
3082e8ab 785
ebdf1dbb 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));
a4c7ef28 798 Float_t phiRel3 = trk->Phi() - fEventPlaneAngle3;
799 if (phiRel3 < 0.)
800 phiRel3 += 2. * TMath::Pi();
ebdf1dbb 801
23536a91 802 if (AcceptAssoc(trk)) {
803 assArray[kAssHad].Add(trk);
804 FillH1(kHistEtaPhiAssHad, trk->Phi(), trk->Eta());
ebdf1dbb 805 FillH2(kHistPhiAssHadEvPlane, phiRel, fCentrality);
a4c7ef28 806 FillH2(kHistPhiAssHadEvPlane3, phiRel3, fCentrality);
3082e8ab 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));
579a71d2 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 }
3082e8ab 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));
579a71d2 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 }
3082e8ab 858 }
23536a91 859 if (IsProton(trk)) {
860 assArray[kAssProt].Add(trk);
861 FillH1(kHistEtaPhiAssProt, trk->Phi(), trk->Eta());
ebdf1dbb 862 FillH2(kHistPhiAssProtEvPlane, phiRel, fCentrality);
a4c7ef28 863 FillH2(kHistPhiAssProtEvPlane3, phiRel3, fCentrality);
23536a91 864 }
3082e8ab 865
866 // template generation
3082e8ab 867 if (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC, trk) == AliPIDResponse::kDetPidOk) {
579a71d2 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
3082e8ab 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);
ebdf1dbb 883
579a71d2 884 if(IsClass(kClCentral))
3082e8ab 885 FillH2(kHistNsigmaTPCe, trk->P(), (rndTPCx-expTPC)/expSigmaTPC, 1., iParticle);
579a71d2 886 if(IsClass(kClSemiCentral))
3082e8ab 887 FillH2(kHistNsigmaTPCeSemi, trk->P(), (rndTPCx-expTPC)/expSigmaTPC, 1., iParticle);
888 }
889 }
890
3082e8ab 891 if (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, trk) == AliPIDResponse::kDetPidOk) {
579a71d2 892 Double_t deltaTOF;
893 fPIDResponse->GetSignalDelta(AliPIDResponse::kTOF, trk, AliPID::kProton, deltaTOF);
894
3082e8ab 895 AliTOFPIDResponse &tofResponse = fPIDResponse->GetTOFResponse();
896
ebdf1dbb 897 Float_t p = trk->P();
50ff8c36 898 if (const AliVVtrack *param = trk->GetInnerParam())
ebdf1dbb 899 p = param->GetP();
3082e8ab 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
ebdf1dbb 910 // Double_t signalSigma = TMath::Sqrt(fTOFsignal.Variance(-2000., 2000.));
3082e8ab 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);
ebdf1dbb 915 // tracking smearing (default parameters, should be overwritten from OADB)
3082e8ab 916 Double_t fPar[] = { 0.008, 0.008, 0.002, 40.0 };
ebdf1dbb 917 if (fParamsTOF)
918 for (Int_t i = 0; i < 4; ++i)
919 fPar[i] = fParamsTOF->GetSigParams(i);
3082e8ab 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));
ebdf1dbb 924 // Double_t cmpSigmaTOFx = fPIDResponse->GetTOFResponse().GetExpectedSigma(p, expTOFx, AliPID::EParticleType(iParticle));
3082e8ab 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));
ebdf1dbb 931 // Double_t tmpSigmaTOFx = TMath::Sqrt(signalSigma*signalSigma +
932 // timezeroSigma*timezeroSigma +
933 // expSigmaTOFx*expSigmaTOFx + fPar[3]*fPar[3]/p/p);
3082e8ab 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);
ebdf1dbb 948 // FillH2(kHistExpSigmaTOFe, p, cmpSigmaTOFx * 1.e-3, 1., iParticle);
949 // FillH2(kHistCmpSigmaTOFe, cmpSigmaTOFx * 1.e-3, tmpSigmaTOFx * 1.e-3, 1., iParticle);
3082e8ab 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);
ebdf1dbb 954 // FillH2(kHistCmpSigmaTOFeSemi, cmpSigmaTOFx * 1.e-3, tmpSigmaTOFx * 1.e-3, 1., iParticle);
955 // FillH2(kHistExpSigmaTOFeSemi, p, cmpSigmaTOFx * 1.e-3, 1., iParticle);
3082e8ab 956 }
957 }
958
959 Double_t rndTOFmismatch = AliTOFPIDResponse::GetMismatchRandomValue(trk->Eta());
960
961 if(IsClass(kClCentral)) {
d57e2346 962 FillH2(kHistNsigmaTOFmismatch, trk->P(), (rndTOFmismatch - expTOF) / expSigmaTOF);
3082e8ab 963 FillH2(kHistDeltaTOF, trk->P(), deltaTOF * 1.e-3); // ps -> ns
964 }
965 if(IsClass(kClSemiCentral)) {
d57e2346 966 FillH2(kHistNsigmaTOFmismatchSemi, trk->P(), (rndTOFmismatch - expTOF) / expSigmaTOF);
3082e8ab 967 FillH2(kHistDeltaTOFSemi, trk->P(), deltaTOF * 1.e-3); // ps -> ns
968 }
969 }
23536a91 970 }
971 }
972
a4c7ef28 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
23536a91 982 // select trigger jet
ebdf1dbb 983 // and remove them from the Q vector
23536a91 984 Int_t nJets = fJetArray ? fJetArray->GetEntries() : 0;
a4c7ef28 985 const TVector2 *qVectorOrig = fEventplane->GetQVector();
986 TVector2 qVector;
987 if (qVectorOrig)
988 qVector = *qVectorOrig;
ebdf1dbb 989
23536a91 990 for (Int_t iJet = 0; iJet < nJets; ++iJet) {
991 AliAODJet *jet = (AliAODJet*) fJetArray->At(iJet);
ebdf1dbb 992
23536a91 993 if (AcceptTrigger(jet)) {
994 trgArray[kTrgJet].Add(jet);
995 FillH1(kHistEtaPhiTrgJet, jet->Phi(), jet->Eta());
ebdf1dbb 996
a4c7ef28 997 if (qVectorOrig) {
ebdf1dbb 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));
a4c7ef28 1005 qVector -= evplaneContrib;
ebdf1dbb 1006 }
1007 }
1008 }
1009 }
1010 }
a4c7ef28 1011 if (qVectorOrig) {
ebdf1dbb 1012 for (Int_t iClass = 0; iClass < kClLast; ++iClass) {
1013 if (IsClass((Class_t) iClass)) {
a4c7ef28 1014 FillH3(kHistEvPlaneCorrNoTrgJets, qVectorOrig->Phi()/2., qVector.Phi()/2., iClass);
ebdf1dbb 1015 if (trgArray[kTrgJet].GetEntriesFast() > 0)
a4c7ef28 1016 FillH3(kHistEvPlaneCorrNoTrgJetsTrgd, qVectorOrig->Phi()/2., qVector.Phi()/2., iClass);
ebdf1dbb 1017 }
23536a91 1018 }
1019 }
1020
ebdf1dbb 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
23536a91 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) {
ebdf1dbb 1034 for (Int_t iAss = 0; iAss <= kAssProt; ++iAss) {
23536a91 1035 // same event
1036 Correlate((Trg_t) iTrg, (Ass_t) iAss, (Class_t) iClass, kEvSame, &trgArray[iTrg], &assArray[iAss]);
1037
1038 // mixed event
2751fb9a 1039 AliEventPool *pool = GetPool((Ass_t) iAss);
23536a91 1040 if (pool && pool->IsReady()) {
2751fb9a 1041 AliDebug(1, Form("----- using pool: %i %i %i -----", iClass, iTrg, iAss));
1042 const Int_t nEvents = pool->GetCurrentNEvents();
23536a91 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 }
23536a91 1048 }
2751fb9a 1049 }
23536a91 1050
2751fb9a 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();
23536a91 1062 }
1063 }
2751fb9a 1064 // }
1065
ebdf1dbb 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]);
23536a91 1071 }
1072 }
ebdf1dbb 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();
23536a91 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
3082e8ab 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
23536a91 1101Bool_t AliAnalysisTaskJetProtonCorr::DetectTriggers()
1102{
1103 fTriggerMask = 0;
1104
ebdf1dbb 1105 AliVEvent::EOfflineTriggerTypes physSel =
1106 (AliVEvent::EOfflineTriggerTypes) ((AliInputEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
579a71d2 1107 // TString trgClasses = InputEvent()->GetFiredTriggerClasses();
23536a91 1108
1109 // physics selection
579a71d2 1110 if (physSel & (AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral))
23536a91 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
ebdf1dbb 1133 // check for run change
1134 if (fRunNumber != InputEvent()->GetRunNumber()) {
1135 fRunNumber = InputEvent()->GetRunNumber();
1136 SetParamsTOF();
1137 }
1138
23536a91 1139 // retrieve z-vertex position
ebdf1dbb 1140 fZvtx = 100.;
23536a91 1141 const AliVVertex *vtx = InputEvent()->GetPrimaryVertex();
ebdf1dbb 1142 if (vtx) {
1143 FillH1(kHistStat, kStatVtx);
1144 FillH1(kHistVertexNctb, vtx->GetNContributors());
1145 if (vtx->GetNContributors() >= 3.)
1146 fZvtx = vtx->GetZ();
1147 }
23536a91 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 }
3082e8ab 1162 else
1163 eventGood = kFALSE;
23536a91 1164
1165 // retrieve event plane
ebdf1dbb 1166 fEventplane = InputEvent()->GetEventplane();
1167 if (fEventplane) {
a4c7ef28 1168 fEventPlaneAngle3 = fEventplane->GetEventplane("V0", InputEvent(), 3);
ebdf1dbb 1169
a4c7ef28 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();
ebdf1dbb 1189 if (fEventPlaneAngleCheck < 0)
1190 fEventPlaneAngleCheck += TMath::Pi();
a4c7ef28 1191
1192 FillH1(kHistStat, kStatEvPlane);
23536a91 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) {
a4c7ef28 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 }
23536a91 1223 }
1224 else if (fAODEvent) {
a4c7ef28 1225 // associate candidates
1226 fPrimTrackArrayAss = new TObjArray();
1227 const Int_t nTracksAODAss = fAODEvent->GetNumberOfTracks();
1228 for (Int_t iTrack = 0; iTrack < nTracksAODAss; ++iTrack) {
23536a91 1229 AliAODTrack *trk = fAODEvent->GetTrack(iTrack);
687dfd6b 1230 if (trk->TestFilterBit(fAssFilterMask))
a4c7ef28 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);
23536a91 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) {
2751fb9a 1250 // still try the output event
1251 if (AODEvent())
1252 fJetArray = dynamic_cast<TClonesArray*> (AODEvent()->FindListObject(fJetBranchName));
1253 if (!fJetArray) {
d57e2346 1254 if (fErrorMsg > 0) {
1255 AliError(Form("no jet branch \"%s\" found, in the AODs are:", fJetBranchName));
1256 --fErrorMsg;
1257 }
2751fb9a 1258 if (fDebug > 0) {
1259 fAODEvent->GetList()->Print();
1260 if (AODEvent())
1261 AODEvent()->GetList()->Print();
1262 }
1263 }
23536a91 1264 }
1265 }
1266
1267 // retrieve event pool for the event category
1268 if (eventGood) {
2751fb9a 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;
23536a91 1274 }
1275 }
1276
1277 return eventGood;
1278}
1279
d57e2346 1280Bool_t AliAnalysisTaskJetProtonCorr::AcceptTrigger(const AliVTrack *trg)
23536a91 1281{
0ba09863 1282 if (TMath::Abs(trg->Eta()) > fHadEtaMax)
1283 return kFALSE;
1284
ebdf1dbb 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();
a4c7ef28 1301 Float_t phiRel3 = trg->Phi() - fEventPlaneAngle3;
1302 if (phiRel3 < 0.)
1303 phiRel3 += 2. * TMath::Pi();
ebdf1dbb 1304 FillH2(kHistPhiTrgHadEvPlane, phiRel, fCentrality);
a4c7ef28 1305 FillH2(kHistPhiTrgHadEvPlane3, phiRel3, fCentrality);
23536a91 1306 }
1307
ebdf1dbb 1308 return (acceptPt && acceptOrientation);
23536a91 1309}
1310
d57e2346 1311AliVTrack* AliAnalysisTaskJetProtonCorr::GetLeadingTrack(const AliAODJet *jet) const
23536a91 1312{
3082e8ab 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;
23536a91 1325 }
1326 }
1327
3082e8ab 1328 // retrieve the leading track
1329 return (AliAODTrack*) jet->GetRefTracks()->At(iLeadingTrack);
23536a91 1330}
1331
d57e2346 1332Bool_t AliAnalysisTaskJetProtonCorr::AcceptTrigger(const AliAODJet *trg)
23536a91 1333{
3082e8ab 1334 // restrict eta
d57e2346 1335 if (TMath::Abs(trg->Eta()) > fTrgJetEtaMax)
1336 return kFALSE;
1337
1338 // reject too small jets
1339 if (trg->EffectiveAreaCharged() < fTrgJetAreaMin)
3082e8ab 1340 return kFALSE;
1341
1342 // require hard leading track
ebdf1dbb 1343 // (leading track biased jet sample)
d57e2346 1344 // but reject jets with too high pt constituents
1345 const Float_t ptLeadTrack = GetLeadingTrack(trg)->Pt();
1346 if ((ptLeadTrack < fTrgJetLeadTrkPtMin) ||
1347 (ptLeadTrack > fTrgJetLeadTrkPtMax))
3082e8ab 1348 return kFALSE;
1349
ebdf1dbb 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;
d57e2346 1355
ebdf1dbb 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();
a4c7ef28 1369 Float_t phiRel3 = trg->Phi() - fEventPlaneAngle3;
1370 if (phiRel3 < 0.)
1371 phiRel3 += 2. * TMath::Pi();
ebdf1dbb 1372 FillH2(kHistPhiTrgJetEvPlane, phiRel, fCentrality);
a4c7ef28 1373 FillH2(kHistPhiTrgJetEvPlane3, phiRel3, fCentrality);
3082e8ab 1374 }
1375
ebdf1dbb 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 }
3082e8ab 1383
ebdf1dbb 1384 return (acceptPt && acceptOrientation);
3082e8ab 1385}
1386
d57e2346 1387Bool_t AliAnalysisTaskJetProtonCorr::AcceptAssoc(const AliVTrack *trk) const
3082e8ab 1388{
687dfd6b 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
3082e8ab 1399 return kTRUE;
687dfd6b 1400 else
1401 return kFALSE;
23536a91 1402}
1403
d57e2346 1404Bool_t AliAnalysisTaskJetProtonCorr::IsProton(const AliVTrack *trk) const
23536a91 1405{
687dfd6b 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
23536a91 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
d57e2346 1421Bool_t AliAnalysisTaskJetProtonCorr::AcceptAngleToEvPlane(Float_t phi, Float_t psi) const
23536a91 1422{
3082e8ab 1423 Float_t deltaPhi = phi - psi;
a4c7ef28 1424 while (deltaPhi < 0.)
1425 deltaPhi += 2 * TMath::Pi();
23536a91 1426
3082e8ab 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);
23536a91 1430
3082e8ab 1431 if (TMath::Abs(deltaPhi) < fTrgAngleToEvPlane)
1432 return kTRUE;
1433 else
1434 return kFALSE;
1435}
23536a91 1436
3082e8ab 1437Float_t AliAnalysisTaskJetProtonCorr::GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1,
1438 Float_t phi2, Float_t pt2, Float_t charge2,
d57e2346 1439 Float_t radius, Float_t bSign) const
3082e8ab 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
d57e2346 1462Bool_t AliAnalysisTaskJetProtonCorr::AcceptTwoTracks(const AliVParticle *trgPart, const AliVParticle *assPart) const
3082e8ab 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
3082e8ab 1474 Float_t deta = trgPart->Eta() - assPart->Eta();
1475
2751fb9a 1476 Float_t bSign = (InputEvent()->GetMagneticField() > 0) ? 1 : -1;
1477
3082e8ab 1478 // optimization
2751fb9a 1479 if (TMath::Abs(deta) < fCutsTwoTrackEff) {
3082e8ab 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) {
2751fb9a 1492 Float_t dphistarabs = TMath::Abs(GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign));
1493 if (dphistarabs < dphistarminabs)
3082e8ab 1494 dphistarminabs = dphistarabs;
3082e8ab 1495 }
1496
2751fb9a 1497 if (dphistarminabs < fCutsTwoTrackEff)
3082e8ab 1498 return kFALSE;
23536a91 1499 }
1500 }
1501
1502 return kTRUE;
1503}
1504
3082e8ab 1505Bool_t AliAnalysisTaskJetProtonCorr::Correlate(CorrType_t corr, Class_t cl, Ev_t ev,
23536a91 1506 TCollection *trgArray, TCollection *assArray, Float_t weight)
1507{
23536a91 1508 AliHistCorr *histCorr = GetHistCorr(corr, cl, ev);
ebdf1dbb 1509 if (!histCorr) {
1510 AliError(Form("no correlation histograms for corr %i, cl %i, ev %i",
1511 corr, cl, ev));
1512 return kFALSE;
1513 }
23536a91 1514
ebdf1dbb 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 }
23536a91 1526
ebdf1dbb 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 }
23536a91 1559 }
1560 }
1561
1562 return kTRUE;
1563}
1564
3082e8ab 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{
3fa79056 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
3082e8ab 1574 CorrType_t corr = (CorrType_t) (2 * trg + ass);
1575
1576 return Correlate(corr, cl, ev, trgArray, assArray, weight);
1577}
1578
ebdf1dbb 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
687dfd6b 1589 const Int_t nPart = gRandom->Poisson(fToyMeanNoPart);
ebdf1dbb 1590
687dfd6b 1591 Float_t trgJetEta = gRandom->Uniform(-fTrgJetEtaMax, fTrgJetEtaMax);
ebdf1dbb 1592 Float_t trgJetPhi = 0.;
0ba09863 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()));
ebdf1dbb 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
687dfd6b 1617 Float_t trgJetEtaAway = gRandom->Uniform(-fHadEtaMax, fHadEtaMax);
1618 Float_t trgJetPhiAway = trgJetPhi + TMath::Pi() + gRandom->Gaus(0., fToySmearPhi);
ebdf1dbb 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;
687dfd6b 1624 const Float_t r = fToyRadius;
ebdf1dbb 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;
3fa79056 1633 if (TMath::Abs(eta) > fHadEtaMax) {
ebdf1dbb 1634 delete assPart;
1635 continue;
1636 }
687dfd6b 1637 while (phi < 0.)
ebdf1dbb 1638 phi += 2. * TMath::Pi();
687dfd6b 1639 while (phi > 2. * TMath::Pi())
ebdf1dbb 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
687dfd6b 1648 Float_t trgHadEta = gRandom->Uniform(-fHadEtaMax, fHadEtaMax);
ebdf1dbb 1649 Float_t trgHadPhi = 0.;
0ba09863 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()));
ebdf1dbb 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
687dfd6b 1674 Float_t trgHadEtaAway = gRandom->Uniform(-fHadEtaMax, fHadEtaMax);
1675 Float_t trgHadPhiAway = trgHadPhi + TMath::Pi() + gRandom->Gaus(0., fToySmearPhi);
ebdf1dbb 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;
687dfd6b 1681 const Float_t r = fToyRadius;
ebdf1dbb 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;
3fa79056 1690 if (TMath::Abs(eta) > fHadEtaMax) {
ebdf1dbb 1691 delete assPart;
1692 continue;
1693 }
687dfd6b 1694 while (phi < 0.)
ebdf1dbb 1695 phi += 2. * TMath::Pi();
687dfd6b 1696 while (phi > 2. * TMath::Pi())
ebdf1dbb 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
23536a91 1707Bool_t AliAnalysisTaskJetProtonCorr::CleanUpEvent()
1708{
1709 if (fAODEvent) {
a4c7ef28 1710 delete fPrimTrackArrayAss;
1711 fPrimTrackArrayAss = 0x0;
1712 delete fPrimTrackArrayTrg;
1713 fPrimTrackArrayTrg = 0x0;
23536a91 1714 }
1715
a4c7ef28 1716 fPrimConstrainedTrackArray->Delete();
1717
23536a91 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++) {
23536a91 1728 // WARNING: TObject::Clone() is very!!! expensive, unusable
1729 // tracksClone->Add(particle->Clone());
1730
0ba09863 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));
23536a91 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),
ebdf1dbb 1749 fHistCorrEtaPhi(0x0),
1750 fHistCorrAvgEtaPhi(0x0),
1751 fHistCorrTrgEtaPhi(0x0),
09fd5f9f 1752 fHistCorrAssEtaPhi(0x0),
1753 fHistDphiLo(-TMath::Pi() / 2.),
1754 fHistDphiNbins(120),
1755 fHistDetaNbins(120)
23536a91 1756{
1757 // ctor
1758
1759 fHistStat = new TH1F(Form("%s_stat", name.Data()), "statistics",
1760 1, .5, 1.5);
1761
09fd5f9f 1762 fHistCorrPhi = new TH1F(Form("%s_phi", name.Data()), ";#Delta#phi",
1763 fHistDphiNbins, fHistDphiLo, 2.*TMath::Pi() + fHistDphiLo);
2751fb9a 1764 fHistCorrPhi->Sumw2();
0ba09863 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();
23536a91 1769 fHistCorrEtaPhi = new TH2F(Form("%s_etaphi", name.Data()), ";#Delta#phi;#Delta#eta",
09fd5f9f 1770 fHistDphiNbins, fHistDphiLo, 2*TMath::Pi() + fHistDphiLo,
1771 fHistDetaNbins, -2., 2.);
2751fb9a 1772 fHistCorrEtaPhi->Sumw2();
0ba09863 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();
23536a91 1785
1786 fOutputList->Add(fHistStat);
1787 fOutputList->Add(fHistCorrPhi);
0ba09863 1788 // fOutputList->Add(fHistCorrPhi2);
23536a91 1789 fOutputList->Add(fHistCorrEtaPhi);
0ba09863 1790 // fOutputList->Add(fHistCorrAvgEtaPhi);
1791 // fOutputList->Add(fHistCorrTrgEtaPhi);
1792 // fOutputList->Add(fHistCorrAssEtaPhi);
23536a91 1793}
1794
1795AliAnalysisTaskJetProtonCorr::AliHistCorr::~AliHistCorr()
1796{
1797 // dtor
1798}
1799
1800void AliAnalysisTaskJetProtonCorr::AliHistCorr::Fill(AliVParticle *trgPart, AliVParticle *assPart, Float_t weight)
1801{
2751fb9a 1802 // fill correlation histograms from trigger particle and associate particle
1803
23536a91 1804 Float_t deltaEta = assPart->Eta() - trgPart->Eta();
ebdf1dbb 1805 Float_t avgEta = (assPart->Eta() + trgPart->Eta()) / 2.;
23536a91 1806 Float_t deltaPhi = assPart->Phi() - trgPart->Phi();
09fd5f9f 1807 while (deltaPhi > (2.*TMath::Pi() + fHistDphiLo))
23536a91 1808 deltaPhi -= 2. * TMath::Pi();
09fd5f9f 1809 while (deltaPhi < fHistDphiLo)
23536a91 1810 deltaPhi += 2. * TMath::Pi();
23536a91 1811
2751fb9a 1812 fHistCorrPhi->Fill(deltaPhi, weight);
0ba09863 1813 if (fHistCorrPhi2)
1814 fHistCorrPhi2->Fill(trgPart->Phi(), assPart->Phi(), weight);
23536a91 1815 fHistCorrEtaPhi->Fill(deltaPhi, deltaEta, weight);
0ba09863 1816 if (fHistCorrAvgEtaPhi)
1817 fHistCorrAvgEtaPhi->Fill(deltaPhi, avgEta, weight);
ebdf1dbb 1818}
1819
1820void AliAnalysisTaskJetProtonCorr::AliHistCorr::Fill(TLorentzVector *trgPart, AliVParticle *assPart, Float_t weight)
1821{
2751fb9a 1822 // fill correlation histograms from trigger direction and associate particle
1823
ebdf1dbb 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();
09fd5f9f 1827 while (deltaPhi > (2.*TMath::Pi() + fHistDphiLo))
ebdf1dbb 1828 deltaPhi -= 2. * TMath::Pi();
09fd5f9f 1829 while (deltaPhi < fHistDphiLo)
ebdf1dbb 1830 deltaPhi += 2. * TMath::Pi();
ebdf1dbb 1831
2751fb9a 1832 fHistCorrPhi->Fill(deltaPhi, weight);
0ba09863 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 }
ebdf1dbb 1839 fHistCorrEtaPhi->Fill(deltaPhi, deltaEta, weight);
0ba09863 1840 if (fHistCorrAvgEtaPhi)
1841 fHistCorrAvgEtaPhi->Fill(deltaPhi, avgEta, weight);
ebdf1dbb 1842}
1843
1844void AliAnalysisTaskJetProtonCorr::AliHistCorr::Fill(TLorentzVector *trgPart, TLorentzVector *assPart, Float_t weight)
1845{
2751fb9a 1846 // fill correlation histograms from trigger direction and associate direction
1847
ebdf1dbb 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();
09fd5f9f 1851 if (deltaPhi > (2.*TMath::Pi() + fHistDphiLo))
ebdf1dbb 1852 deltaPhi -= 2. * TMath::Pi();
09fd5f9f 1853 else if (deltaPhi < fHistDphiLo)
ebdf1dbb 1854 deltaPhi += 2. * TMath::Pi();
ebdf1dbb 1855
2751fb9a 1856 fHistCorrPhi->Fill(deltaPhi, weight);
0ba09863 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 }
ebdf1dbb 1866 fHistCorrEtaPhi->Fill(deltaPhi, deltaEta, weight);
0ba09863 1867 if (fHistCorrAvgEtaPhi)
1868 fHistCorrAvgEtaPhi->Fill(deltaPhi, avgEta, weight);
23536a91 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}