]> 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
- V0base: Implemented different vertex cut: also accept SPD vertices in ESDs. For...
[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) {
23536a91 1248 AliAODTrack *trk = fAODEvent->GetTrack(iTrack);
687dfd6b 1249 if (trk->TestFilterBit(fAssFilterMask))
a4c7ef28 1250 fPrimTrackArrayAss->Add(trk);
1251 }
1252
1253 // trigger candidates
1254 fPrimTrackArrayTrg = new TObjArray();
1255 const Int_t nTracksAODTrg = fAODEvent->GetNumberOfTracks();
1256 for (Int_t iTrack = 0; iTrack < nTracksAODTrg; ++iTrack) {
1257 AliAODTrack *trk = fAODEvent->GetTrack(iTrack);
1258 if (trk->IsHybridGlobalConstrainedGlobal())
1259 fPrimTrackArrayTrg->Add(trk);
23536a91 1260 }
1261 }
1262 else
1263 eventGood = kFALSE;
1264
1265 // retrieve jet array
1266 if (fAODEvent) {
1267 fJetArray = dynamic_cast<TClonesArray*> (fAODEvent->FindListObject(fJetBranchName));
1268 if (!fJetArray) {
2751fb9a 1269 // still try the output event
1270 if (AODEvent())
1271 fJetArray = dynamic_cast<TClonesArray*> (AODEvent()->FindListObject(fJetBranchName));
1272 if (!fJetArray) {
d57e2346 1273 if (fErrorMsg > 0) {
1274 AliError(Form("no jet branch \"%s\" found, in the AODs are:", fJetBranchName));
1275 --fErrorMsg;
1276 }
2751fb9a 1277 if (fDebug > 0) {
1278 fAODEvent->GetList()->Print();
1279 if (AODEvent())
1280 AODEvent()->GetList()->Print();
1281 }
1282 }
23536a91 1283 }
1284 }
1285
1286 // retrieve event pool for the event category
1287 if (eventGood) {
2751fb9a 1288 for (Int_t iAss = 0; iAss <= kAssProt; ++iAss) {
1289 AliEventPoolManager *mgr = GetPoolMgr((Ass_t) iAss);
1290 GetPool((Ass_t) iAss) =
1291 // mgr ? mgr->GetEventPool(fCentrality, fZvtx, fEventPlaneAngle) : 0x0;
1292 mgr ? mgr->GetEventPool(fCentrality, fZvtx) : 0x0;
23536a91 1293 }
1294 }
1295
1296 return eventGood;
1297}
1298
af846fc0 1299Float_t AliAnalysisTaskJetProtonCorr::GetPhiRel2(AliVParticle *part) const
1300{
1301 // calculate the angle to the event plane
1302 // after removing the particle's own contribution to the Q vector
1303
1304 AliVTrack *track = dynamic_cast<AliVTrack*> (part);
1305 AliAODJet *jet = dynamic_cast<AliAODJet*> (part);
1306
1307 Float_t phiRel = -1.;
1308
1309 if (!fEventplane)
1310 return phiRel;
1311
1312 const TVector2 *qVectorOrig = fEventplane->GetQVector();
1313 if (!qVectorOrig)
1314 return phiRel;
1315
1316 TVector2 qVector = *qVectorOrig;
1317
50e3095f 1318 // ??? protect against negative ID
1319 // but should be handled properly by event plane
1320 if (track && (track->GetID() > -1)) {
af846fc0 1321 TVector2 evplaneContrib(fEventplane->GetQContributionX(track),
1322 fEventplane->GetQContributionY(track));
1323 qVector -= evplaneContrib;
1324 }
1325 else if (jet) {
1326 const Int_t nRefTracks = jet->GetRefTracks()->GetEntriesFast();
1327 for (Int_t iTrack = 0; iTrack < nRefTracks; ++iTrack) {
1328 AliVTrack *contrib = (AliVTrack*) jet->GetRefTracks()->At(iTrack);
1329
50e3095f 1330 if (contrib && (contrib->GetID() > -1)) {
1331 TVector2 evplaneContrib(fEventplane->GetQContributionX(contrib),
1332 fEventplane->GetQContributionY(contrib));
1333 qVector -= evplaneContrib;
af846fc0 1334 }
1335 }
1336 }
1337
1338 // assign phiRel and map to [0, 2 pi)
1339 phiRel = part->Phi() - (qVector.Phi() / 2.);
1340 // if (gRandom->Rndm() > .5)
1341 // phiRel -= TMath::Pi();
1342 while (phiRel < 0.)
1343 phiRel += TMath::TwoPi();
1344 while (phiRel > TMath::TwoPi())
1345 phiRel -= TMath::TwoPi();
1346
1347 return phiRel;
1348}
1349
1350Bool_t AliAnalysisTaskJetProtonCorr::AcceptTrigger(AliVTrack *trg)
23536a91 1351{
0ba09863 1352 if (TMath::Abs(trg->Eta()) > fHadEtaMax)
1353 return kFALSE;
1354
ebdf1dbb 1355 // check pt interval
1356 Bool_t acceptPt =
1357 (trg->Pt() >= fTrgPartPtMin) &&
1358 (trg->Pt() <= fTrgPartPtMax);
1359
1360 // for semi-central events check phi w.r.t. event plane
1361 Bool_t acceptOrientation = IsSemiCentral() ?
1362 AcceptAngleToEvPlane(trg->Phi(), GetEventPlaneAngle()) :
1363 kTRUE;
1364
1365 if (acceptPt) {
1366 Float_t phiRel = trg->Phi() - fEventPlaneAngle;
af846fc0 1367 // if (gRandom->Rndm() > .5)
1368 // phiRel -= TMath::Pi();
ebdf1dbb 1369 if (phiRel < 0.)
1370 phiRel += 2. * TMath::Pi();
a4c7ef28 1371 Float_t phiRel3 = trg->Phi() - fEventPlaneAngle3;
1372 if (phiRel3 < 0.)
1373 phiRel3 += 2. * TMath::Pi();
af846fc0 1374 FillH2(kHistPhiTrgHadEvPlane, GetPhiRel2(trg), fCentrality);
a4c7ef28 1375 FillH2(kHistPhiTrgHadEvPlane3, phiRel3, fCentrality);
23536a91 1376 }
1377
ebdf1dbb 1378 return (acceptPt && acceptOrientation);
23536a91 1379}
1380
d57e2346 1381AliVTrack* AliAnalysisTaskJetProtonCorr::GetLeadingTrack(const AliAODJet *jet) const
23536a91 1382{
3082e8ab 1383 // return leading track within a jet
1384
1385 // check contributing tracks
1386 Int_t nJetTracks = jet->GetRefTracks()->GetEntriesFast();
1387 Int_t iLeadingTrack = -1;
1388 Float_t ptLeadingTrack = 0.;
1389 for (Int_t iTrack=0; iTrack < nJetTracks; ++iTrack) {
1390 AliAODTrack *track = (AliAODTrack*) jet->GetRefTracks()->At(iTrack);
1391 // find the leading track
1392 if (track->Pt() > ptLeadingTrack) {
1393 ptLeadingTrack = track->Pt();
1394 iLeadingTrack = iTrack;
23536a91 1395 }
1396 }
1397
3082e8ab 1398 // retrieve the leading track
1399 return (AliAODTrack*) jet->GetRefTracks()->At(iLeadingTrack);
23536a91 1400}
1401
af846fc0 1402Bool_t AliAnalysisTaskJetProtonCorr::AcceptTrigger(AliAODJet *trg)
23536a91 1403{
3082e8ab 1404 // restrict eta
d57e2346 1405 if (TMath::Abs(trg->Eta()) > fTrgJetEtaMax)
1406 return kFALSE;
1407
1408 // reject too small jets
1409 if (trg->EffectiveAreaCharged() < fTrgJetAreaMin)
3082e8ab 1410 return kFALSE;
1411
1412 // require hard leading track
ebdf1dbb 1413 // (leading track biased jet sample)
d57e2346 1414 // but reject jets with too high pt constituents
1415 const Float_t ptLeadTrack = GetLeadingTrack(trg)->Pt();
1416 if ((ptLeadTrack < fTrgJetLeadTrkPtMin) ||
1417 (ptLeadTrack > fTrgJetLeadTrkPtMax))
3082e8ab 1418 return kFALSE;
1419
ebdf1dbb 1420 // check for jet orientation w.r.t. event plane
1421 // (constrained only for semi-central events)
1422 Bool_t acceptOrientation = IsSemiCentral() ?
1423 AcceptAngleToEvPlane(trg->Phi(), GetEventPlaneAngle()) :
1424 kTRUE;
d57e2346 1425
ebdf1dbb 1426 // check for jet pt
1427 Bool_t acceptPt =
1428 (trg->Pt() >= fTrgJetPtMin) &&
1429 (trg->Pt() <= fTrgJetPtMax);
1430
1431 if (acceptPt) {
1432 // store the azimuthal distribution relative to the event plane
1433 // for the jets in the pt range of interest
1434 Float_t phiRel = trg->Phi() - fEventPlaneAngle;
af846fc0 1435 // if (gRandom->Rndm() > .5)
1436 // phiRel -= TMath::Pi();
ebdf1dbb 1437 if (phiRel < 0.)
1438 phiRel += 2. * TMath::Pi();
a4c7ef28 1439 Float_t phiRel3 = trg->Phi() - fEventPlaneAngle3;
1440 if (phiRel3 < 0.)
1441 phiRel3 += 2. * TMath::Pi();
af846fc0 1442 FillH2(kHistPhiTrgJetEvPlane, GetPhiRel2(trg), fCentrality);
a4c7ef28 1443 FillH2(kHistPhiTrgJetEvPlane3, phiRel3, fCentrality);
3082e8ab 1444 }
1445
ebdf1dbb 1446 if (acceptOrientation) {
1447 // spectrum for cross checks
1448 if (IsClass(kClCentral))
1449 FillH1(kHistJetPtCentral, trg->Pt());
1450 if (IsClass(kClSemiCentral))
1451 FillH1(kHistJetPtSemi, trg->Pt());
1452 }
3082e8ab 1453
ebdf1dbb 1454 return (acceptPt && acceptOrientation);
3082e8ab 1455}
1456
d57e2346 1457Bool_t AliAnalysisTaskJetProtonCorr::AcceptAssoc(const AliVTrack *trk) const
3082e8ab 1458{
687dfd6b 1459 if ((trk->Pt() > fAssPartPtMin) && (trk->Pt() < fAssPartPtMax) && (TMath::Abs(trk->Eta()) < fHadEtaMax))
1460 if (fRequirePID) {
1461 if ((fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC, trk) == AliPIDResponse::kDetPidOk) &&
1462 (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, trk) == AliPIDResponse::kDetPidOk) &&
1463 (trk->GetTPCsignalN() >= 60.) && (trk->GetTPCsignal() >= 10))
1464 return kTRUE;
1465 else
1466 return kFALSE;
1467 }
1468 else
3082e8ab 1469 return kTRUE;
687dfd6b 1470 else
1471 return kFALSE;
23536a91 1472}
1473
d57e2346 1474Bool_t AliAnalysisTaskJetProtonCorr::IsProton(const AliVTrack *trk) const
23536a91 1475{
687dfd6b 1476 if ((fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC, trk) != AliPIDResponse::kDetPidOk) ||
1477 (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, trk) != AliPIDResponse::kDetPidOk) ||
1478 (trk->GetTPCsignalN() < 60.) || (trk->GetTPCsignal() < 10))
1479 return kFALSE;
1480
23536a91 1481 Double_t nSigmaProtonTPC = fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton);
1482 Double_t nSigmaProtonTOF = fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton);
1483
1484 if ((TMath::Abs(nSigmaProtonTPC) <= 2.) && (TMath::Abs(nSigmaProtonTOF) <= 2.)) {
1485 return kTRUE;
1486 }
1487
1488 return kFALSE;
1489}
1490
d57e2346 1491Bool_t AliAnalysisTaskJetProtonCorr::AcceptAngleToEvPlane(Float_t phi, Float_t psi) const
23536a91 1492{
3082e8ab 1493 Float_t deltaPhi = phi - psi;
a4c7ef28 1494 while (deltaPhi < 0.)
1495 deltaPhi += 2 * TMath::Pi();
23536a91 1496
3082e8ab 1497 // map to interval [-pi/2, pi/2)
1498 deltaPhi = std::fmod(deltaPhi + TMath::Pi()/2., TMath::Pi()) - TMath::Pi()/2.;
1499 // printf("delta phi = %f with phi = %f, psi = %f\n", deltaPhi, phi, psi);
23536a91 1500
3082e8ab 1501 if (TMath::Abs(deltaPhi) < fTrgAngleToEvPlane)
1502 return kTRUE;
1503 else
1504 return kFALSE;
1505}
23536a91 1506
3082e8ab 1507Float_t AliAnalysisTaskJetProtonCorr::GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1,
1508 Float_t phi2, Float_t pt2, Float_t charge2,
d57e2346 1509 Float_t radius, Float_t bSign) const
3082e8ab 1510{
1511 // calculates dphistar
1512 // from AliUEHistograms
1513
1514 Float_t dphistar = phi1 - phi2
1515 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1)
1516 + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
1517
1518 static const Double_t kPi = TMath::Pi();
1519
1520 // circularity
1521 if (dphistar > kPi)
1522 dphistar = kPi * 2 - dphistar;
1523 if (dphistar < -kPi)
1524 dphistar = -kPi * 2 - dphistar;
1525 if (dphistar > kPi) // might look funny but is needed
1526 dphistar = kPi * 2 - dphistar;
1527
1528 return dphistar;
1529}
1530
1531
d57e2346 1532Bool_t AliAnalysisTaskJetProtonCorr::AcceptTwoTracks(const AliVParticle *trgPart, const AliVParticle *assPart) const
3082e8ab 1533{
1534 // apply two track pair cut
1535
1536 Float_t phi1 = trgPart->Phi();
1537 Float_t pt1 = trgPart->Pt();
1538 Float_t charge1 = trgPart->Charge();
1539
1540 Float_t phi2 = assPart->Phi();
1541 Float_t pt2 = assPart->Pt();
1542 Float_t charge2 = assPart->Charge();
1543
3082e8ab 1544 Float_t deta = trgPart->Eta() - assPart->Eta();
1545
2751fb9a 1546 Float_t bSign = (InputEvent()->GetMagneticField() > 0) ? 1 : -1;
1547
3082e8ab 1548 // optimization
2751fb9a 1549 if (TMath::Abs(deta) < fCutsTwoTrackEff) {
3082e8ab 1550 // check first boundaries to see if is worth to loop and find the minimum
1551 Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 0.8, bSign);
1552 Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 2.5, bSign);
1553
1554 const Float_t kLimit = fCutsTwoTrackEff * 3;
1555
1556 Float_t dphistarminabs = 1e5;
1557 // Float_t dphistarmin = 1e5;
1558 if ((TMath::Abs(dphistar1) < kLimit) ||
1559 (TMath::Abs(dphistar2) < kLimit) ||
1560 (dphistar1 * dphistar2 < 0)) {
1561 for (Double_t rad=0.8; rad<2.51; rad+=0.01) {
2751fb9a 1562 Float_t dphistarabs = TMath::Abs(GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign));
1563 if (dphistarabs < dphistarminabs)
3082e8ab 1564 dphistarminabs = dphistarabs;
3082e8ab 1565 }
1566
2751fb9a 1567 if (dphistarminabs < fCutsTwoTrackEff)
3082e8ab 1568 return kFALSE;
23536a91 1569 }
1570 }
1571
1572 return kTRUE;
1573}
1574
3082e8ab 1575Bool_t AliAnalysisTaskJetProtonCorr::Correlate(CorrType_t corr, Class_t cl, Ev_t ev,
23536a91 1576 TCollection *trgArray, TCollection *assArray, Float_t weight)
1577{
23536a91 1578 AliHistCorr *histCorr = GetHistCorr(corr, cl, ev);
ebdf1dbb 1579 if (!histCorr) {
1580 AliError(Form("no correlation histograms for corr %i, cl %i, ev %i",
1581 corr, cl, ev));
1582 return kFALSE;
1583 }
23536a91 1584
ebdf1dbb 1585 TIter assIter(assArray);
1586 while (TObject *assObj = assIter()) {
1587 if (dynamic_cast<AliVParticle*> (assObj)) {
1588 AliVParticle *assPart = (AliVParticle*) assObj;
af846fc0 1589 histCorr->Ass(assPart->Phi(), assPart->Eta(),
1590 assPart->Pt() > 0. ? (assPart->Charge() / 3.) / assPart->Pt() : 1.e-4,
1591 weight);
ebdf1dbb 1592 }
1593 else if (dynamic_cast<TLorentzVector*> (assObj)) {
1594 TLorentzVector *assVec = (TLorentzVector*) assObj;
50e3095f 1595 histCorr->Ass(TVector2::Phi_0_2pi(assVec->Phi()), assVec->Eta(), 0., weight);
ebdf1dbb 1596 }
1597 }
23536a91 1598
ebdf1dbb 1599 TIter trgIter(trgArray);
1600 while (TObject *trgObj = trgIter()) {
1601 if (AliVParticle *trgPart = dynamic_cast<AliVParticle*> (trgObj)) {
1602 // count the trigger
af846fc0 1603 histCorr->Trigger(trgPart->Phi(), trgPart->Eta(),
1604 trgPart->Pt() > 0. ? (trgPart->Charge() / 3.) / trgPart->Pt() : 1.e-4,
1605 weight);
ebdf1dbb 1606
1607 // loop over associates
1608 assIter.Reset();
1609 while (TObject *assObj = assIter()) {
1610 if (AliVParticle *assPart = dynamic_cast<AliVParticle*> (assObj)) {
1611 if (AcceptTwoTracks(trgPart, assPart))
1612 histCorr->Fill(trgPart, assPart, weight);
1613 }
1614 else if (TLorentzVector *assVec = dynamic_cast<TLorentzVector*> (assObj)) {
1615 AliFatal(Form("got %p, but not implemented", assVec));
1616 }
1617 }
1618 }
1619 else if (TLorentzVector *trgVec = dynamic_cast<TLorentzVector*> (trgObj)) {
1620 // count the trigger
50e3095f 1621 histCorr->Trigger(TVector2::Phi_0_2pi(trgVec->Phi()), trgVec->Eta(),
af846fc0 1622 0.,
1623 weight);
ebdf1dbb 1624
1625 // loop over associates
1626 assIter.Reset();
1627 while (TObject *assObj = assIter()) {
1628 if (AliVParticle *assPart = dynamic_cast<AliVParticle*> (assObj)) {
1629 histCorr->Fill(trgVec, assPart, weight);
1630 }
1631 else if (TLorentzVector *assVec = dynamic_cast<TLorentzVector*> (assObj)) {
1632 histCorr->Fill(trgVec, assVec, weight);
1633 }
1634 }
23536a91 1635 }
1636 }
1637
1638 return kTRUE;
1639}
1640
3082e8ab 1641Bool_t AliAnalysisTaskJetProtonCorr::Correlate(Trg_t trg, Ass_t ass, Class_t cl, Ev_t ev,
1642 TCollection *trgArray, TCollection *assArray, Float_t weight)
1643{
3fa79056 1644 if ((trg < 0) || (trg > kTrgJetRnd))
1645 AliFatal(Form("wrong request for correlation with trigger: %d", trg));
1646
1647 if ((ass < 0) || (ass > kAssProt))
1648 AliFatal(Form("wrong request for correlation with associate: %d", ass));
1649
3082e8ab 1650 CorrType_t corr = (CorrType_t) (2 * trg + ass);
1651
1652 return Correlate(corr, cl, ev, trgArray, assArray, weight);
1653}
1654
ebdf1dbb 1655Bool_t AliAnalysisTaskJetProtonCorr::GenerateRandom(TCollection *trgJetArray,
1656 TCollection *trgHadArray,
1657 TCollection *assHadJetArray,
1658 TCollection *assProtJetArray,
1659 TCollection *assHadHadArray,
1660 TCollection *assProtHadArray,
1661 Float_t pFraction)
1662{
687dfd6b 1663 const Int_t nPart = gRandom->Poisson(fToyMeanNoPart);
ebdf1dbb 1664
af846fc0 1665 // overcompensate modulation for event plane resolution
1666 if (fSplineEventPlaneRes) {
1667 fTrgJetPhiModCent->SetParameter(0, fTrgJetV2Cent / fSplineEventPlaneRes->Eval(GetCentrality()));
1668 fTrgJetPhiModSemi->SetParameter(0, fTrgJetV2Semi / fSplineEventPlaneRes->Eval(GetCentrality()));
1669 fTrgHadPhiModCent->SetParameter(0, fTrgHadV2Cent / fSplineEventPlaneRes->Eval(GetCentrality()));
1670 fTrgHadPhiModSemi->SetParameter(0, fTrgHadV2Semi / fSplineEventPlaneRes->Eval(GetCentrality()));
1671 }
1672 else {
1673 printf("no event plane resolution spline available\n");
1674 fTrgJetPhiModCent->SetParameter(0, fTrgJetV2Cent);
1675 fTrgJetPhiModSemi->SetParameter(0, fTrgJetV2Semi);
1676 fTrgHadPhiModCent->SetParameter(0, fTrgHadV2Cent);
1677 fTrgHadPhiModSemi->SetParameter(0, fTrgHadV2Semi);
1678 }
1679
1680 // generate random direction
687dfd6b 1681 Float_t trgJetEta = gRandom->Uniform(-fTrgJetEtaMax, fTrgJetEtaMax);
ebdf1dbb 1682 Float_t trgJetPhi = 0.;
0ba09863 1683 do {
1684 trgJetPhi = IsClass(kClSemiCentral) ?
1685 fTrgJetPhiModSemi->GetRandom() :
1686 fTrgJetPhiModCent->GetRandom();
1687
1688 trgJetPhi += fEventPlaneAngle;
1689 if (trgJetPhi < 0.)
1690 trgJetPhi += 2. * TMath::Pi();
1691 else if (trgJetPhi > 2. * TMath::Pi())
1692 trgJetPhi -= 2. * TMath::Pi();
1693
1694 Float_t phiRel = trgJetPhi- fEventPlaneAngle;
1695 if (phiRel < 0.)
1696 phiRel += 2. * TMath::Pi();
1697 FillH2(kHistPhiRndTrgJetEvPlane, phiRel, fCentrality);
1698 } while (IsClass(kClSemiCentral) &&
1699 !AcceptAngleToEvPlane(trgJetPhi, GetEventPlaneAngle()));
ebdf1dbb 1700
1701 // generate one trigger particle
1702 TLorentzVector *trgJet = new TLorentzVector();
1703 trgJet->SetPtEtaPhiM(50., trgJetEta, trgJetPhi, .14);
1704 trgJetArray->Add(trgJet);
1705
1706 // generate direction for away side
687dfd6b 1707 Float_t trgJetEtaAway = gRandom->Uniform(-fHadEtaMax, fHadEtaMax);
1708 Float_t trgJetPhiAway = trgJetPhi + TMath::Pi() + gRandom->Gaus(0., fToySmearPhi);
ebdf1dbb 1709
1710 // generate associated particles
1711 // with proton/hadron ratio observed in this event
1712 for (Int_t iPart = 0; iPart < nPart; ++iPart) {
1713 Float_t deltaEta, deltaPhi;
687dfd6b 1714 const Float_t r = fToyRadius;
ebdf1dbb 1715 do {
1716 deltaEta = gRandom->Uniform(-r, r);
1717 deltaPhi = gRandom->Uniform(-r, r);
1718 } while ((deltaEta * deltaEta + deltaPhi * deltaPhi) > (r * r));
1719
1720 TLorentzVector *assPart = new TLorentzVector();
1721 Float_t eta = trgJetEtaAway + deltaEta;
1722 Float_t phi = trgJetPhiAway + deltaPhi;
3fa79056 1723 if (TMath::Abs(eta) > fHadEtaMax) {
ebdf1dbb 1724 delete assPart;
1725 continue;
1726 }
687dfd6b 1727 while (phi < 0.)
ebdf1dbb 1728 phi += 2. * TMath::Pi();
687dfd6b 1729 while (phi > 2. * TMath::Pi())
ebdf1dbb 1730 phi -= 2. * TMath::Pi();
1731 assPart->SetPtEtaPhiM(3., eta, phi, .14);
1732 assHadJetArray->Add(assPart);
1733 if (gRandom->Uniform() < pFraction)
1734 assProtJetArray->Add(assPart);
1735 }
1736
1737 // trigger hadron
687dfd6b 1738 Float_t trgHadEta = gRandom->Uniform(-fHadEtaMax, fHadEtaMax);
ebdf1dbb 1739 Float_t trgHadPhi = 0.;
0ba09863 1740 do {
1741 trgHadPhi = IsClass(kClSemiCentral) ?
1742 fTrgHadPhiModSemi->GetRandom() :
1743 fTrgHadPhiModCent->GetRandom();
1744
1745 trgHadPhi += fEventPlaneAngle;
1746 if (trgHadPhi < 0.)
1747 trgHadPhi += 2. * TMath::Pi();
1748 else if (trgHadPhi > 2. * TMath::Pi())
1749 trgHadPhi -= 2. * TMath::Pi();
1750
1751 Float_t phiRel = trgHadPhi - fEventPlaneAngle;
1752 if (phiRel < 0.)
1753 phiRel += 2. * TMath::Pi();
1754 FillH2(kHistPhiRndTrgHadEvPlane, phiRel, fCentrality);
1755 } while (IsClass(kClSemiCentral) &&
1756 !AcceptAngleToEvPlane(trgHadPhi, GetEventPlaneAngle()));
ebdf1dbb 1757
1758 // generate one trigger particle
1759 TLorentzVector *trgHad = new TLorentzVector();
1760 trgHad->SetPtEtaPhiM(50., trgHadEta, trgHadPhi, .14);
1761 trgHadArray->Add(trgHad);
1762
1763 // generate direction for away side
687dfd6b 1764 Float_t trgHadEtaAway = gRandom->Uniform(-fHadEtaMax, fHadEtaMax);
1765 Float_t trgHadPhiAway = trgHadPhi + TMath::Pi() + gRandom->Gaus(0., fToySmearPhi);
ebdf1dbb 1766
1767 // generate associated particles
1768 // with proton/hadron ratio observed in this event
1769 for (Int_t iPart = 0; iPart < nPart; ++iPart) {
1770 Float_t deltaEta, deltaPhi;
687dfd6b 1771 const Float_t r = fToyRadius;
ebdf1dbb 1772 do {
1773 deltaEta = gRandom->Uniform(-r, r);
1774 deltaPhi = gRandom->Uniform(-r, r);
1775 } while ((deltaEta * deltaEta + deltaPhi * deltaPhi) > (r * r));
1776
1777 TLorentzVector *assPart = new TLorentzVector();
1778 Float_t eta = trgHadEtaAway + deltaEta;
1779 Float_t phi = trgHadPhiAway + deltaPhi;
3fa79056 1780 if (TMath::Abs(eta) > fHadEtaMax) {
ebdf1dbb 1781 delete assPart;
1782 continue;
1783 }
687dfd6b 1784 while (phi < 0.)
ebdf1dbb 1785 phi += 2. * TMath::Pi();
687dfd6b 1786 while (phi > 2. * TMath::Pi())
ebdf1dbb 1787 phi -= 2. * TMath::Pi();
1788 assPart->SetPtEtaPhiM(3., eta, phi, .14);
1789 assHadHadArray->Add(assPart);
1790 if (gRandom->Uniform() < pFraction)
1791 assProtHadArray->Add(assPart);
1792 }
1793
1794 return kTRUE;
1795}
1796
23536a91 1797Bool_t AliAnalysisTaskJetProtonCorr::CleanUpEvent()
1798{
1799 if (fAODEvent) {
a4c7ef28 1800 delete fPrimTrackArrayAss;
1801 fPrimTrackArrayAss = 0x0;
1802 delete fPrimTrackArrayTrg;
1803 fPrimTrackArrayTrg = 0x0;
23536a91 1804 }
1805
a4c7ef28 1806 fPrimConstrainedTrackArray->Delete();
1807
23536a91 1808 return kTRUE;
1809}
1810
1811TObjArray* AliAnalysisTaskJetProtonCorr::CloneTracks(TObjArray* tracks) const
1812{
1813 TObjArray* tracksClone = new TObjArray;
1814 tracksClone->SetOwner(kTRUE);
1815
1816 Int_t nTracks = tracks->GetEntriesFast();
1817 for (Int_t i = 0; i < nTracks; i++) {
23536a91 1818 // WARNING: TObject::Clone() is very!!! expensive, unusable
1819 // tracksClone->Add(particle->Clone());
1820
0ba09863 1821 // if (AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (tracks->At(i)))
1822 // tracksClone->Add(new AliESDtrack(*esdTrack));
1823 // else if (AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*> (tracks->At(i)))
1824 // tracksClone->Add(new AliAODTrack(*aodTrack));
1825
1826 if (const AliVTrack *track = dynamic_cast<const AliVTrack*> (tracks->At(i)))
1827 tracksClone->Add(new AliPartCorr(*track));
23536a91 1828 }
1829
1830 return tracksClone;
1831}
1832
1833AliAnalysisTaskJetProtonCorr::AliHistCorr::AliHistCorr(TString name, TList *outputList) :
1834 TNamed(name, name),
1835 fOutputList(outputList),
1836 fHistStat(0x0),
1837 fHistCorrPhi(0x0),
1838 fHistCorrPhi2(0x0),
ebdf1dbb 1839 fHistCorrEtaPhi(0x0),
1840 fHistCorrAvgEtaPhi(0x0),
1841 fHistCorrTrgEtaPhi(0x0),
09fd5f9f 1842 fHistCorrAssEtaPhi(0x0),
af846fc0 1843 fHistCorrTrgEtaPhiQpt(0x0),
1844 fHistCorrAssEtaPhiQpt(0x0),
09fd5f9f 1845 fHistDphiLo(-TMath::Pi() / 2.),
1846 fHistDphiNbins(120),
1847 fHistDetaNbins(120)
23536a91 1848{
1849 // ctor
1850
1851 fHistStat = new TH1F(Form("%s_stat", name.Data()), "statistics",
1852 1, .5, 1.5);
1853
af846fc0 1854 fHistCorrPhi = new TH1F(Form("%s_phi", name.Data()), ";#Delta#varphi",
09fd5f9f 1855 fHistDphiNbins, fHistDphiLo, 2.*TMath::Pi() + fHistDphiLo);
2751fb9a 1856 fHistCorrPhi->Sumw2();
af846fc0 1857 fHistCorrPhi2 = new TH2F(Form("%s_phi2", name.Data()), ";#varphi_{trg};#varphi_{ass}",
1858 120, 0.*TMath::Pi(), 2.*TMath::Pi(),
1859 120, 0.*TMath::Pi(), 2.*TMath::Pi());
1860 fHistCorrPhi2->Sumw2();
1861 fHistCorrEtaPhi = new TH2F(Form("%s_etaphi", name.Data()), ";#Delta#varphi;#Delta#eta",
09fd5f9f 1862 fHistDphiNbins, fHistDphiLo, 2*TMath::Pi() + fHistDphiLo,
1863 fHistDetaNbins, -2., 2.);
2751fb9a 1864 fHistCorrEtaPhi->Sumw2();
af846fc0 1865 fHistCorrAvgEtaPhi = new TH2F(Form("%s_avgetaphi", name.Data()), ";#Delta#varphi;avg #eta",
1866 fHistDphiNbins, fHistDphiLo, 2*TMath::Pi() + fHistDphiLo,
1867 fHistDetaNbins, -2., 2.);
1868 fHistCorrAvgEtaPhi->Sumw2();
1869 fHistCorrTrgEtaPhi = new TH2F(Form("%s_trg_etaphi", name.Data()), ";#varphi;#eta",
1870 120, 0., 2*TMath::Pi(),
1871 120, -1., 1.);
1872 fHistCorrTrgEtaPhi->Sumw2();
1873 fHistCorrAssEtaPhi = new TH2F(Form("%s_ass_etaphi", name.Data()), ";#varphi;#eta",
1874 120, 0., 2*TMath::Pi(),
1875 120, -1., 1.);
1876 fHistCorrAssEtaPhi->Sumw2();
1877
1878 // fHistCorrTrgEtaPhiQpt = new TH3F(Form("%s_trg_etaphi", name.Data()), ";#varphi;#eta;Q/p_{T}",
1879 // 120, 0., 2*TMath::Pi(),
1880 // 120, -1., 1.,
1881 // 120, -1., 1.);
1882 // fHistCorrTrgEtaPhiQpt->Sumw2();
1883 // fHistCorrAssEtaPhiQpt = new TH3F(Form("%s_ass_etaphi", name.Data()), ";#varphi;#etaQ/p_{T}",
1884 // 120, 0., 2*TMath::Pi(),
1885 // 120, -1., 1.,
1886 // 120, -1., 1.);
1887 // fHistCorrAssEtaPhiQpt->Sumw2();
23536a91 1888
1889 fOutputList->Add(fHistStat);
1890 fOutputList->Add(fHistCorrPhi);
af846fc0 1891 fOutputList->Add(fHistCorrPhi2);
23536a91 1892 fOutputList->Add(fHistCorrEtaPhi);
af846fc0 1893 fOutputList->Add(fHistCorrAvgEtaPhi);
1894 fOutputList->Add(fHistCorrTrgEtaPhi);
1895 fOutputList->Add(fHistCorrAssEtaPhi);
1896 // fOutputList->Add(fHistCorrTrgEtaPhiQpt);
1897 // fOutputList->Add(fHistCorrAssEtaPhiQpt);
23536a91 1898}
1899
1900AliAnalysisTaskJetProtonCorr::AliHistCorr::~AliHistCorr()
1901{
1902 // dtor
1903}
1904
1905void AliAnalysisTaskJetProtonCorr::AliHistCorr::Fill(AliVParticle *trgPart, AliVParticle *assPart, Float_t weight)
1906{
2751fb9a 1907 // fill correlation histograms from trigger particle and associate particle
1908
23536a91 1909 Float_t deltaEta = assPart->Eta() - trgPart->Eta();
ebdf1dbb 1910 Float_t avgEta = (assPart->Eta() + trgPart->Eta()) / 2.;
23536a91 1911 Float_t deltaPhi = assPart->Phi() - trgPart->Phi();
09fd5f9f 1912 while (deltaPhi > (2.*TMath::Pi() + fHistDphiLo))
23536a91 1913 deltaPhi -= 2. * TMath::Pi();
09fd5f9f 1914 while (deltaPhi < fHistDphiLo)
23536a91 1915 deltaPhi += 2. * TMath::Pi();
23536a91 1916
2751fb9a 1917 fHistCorrPhi->Fill(deltaPhi, weight);
0ba09863 1918 if (fHistCorrPhi2)
1919 fHistCorrPhi2->Fill(trgPart->Phi(), assPart->Phi(), weight);
23536a91 1920 fHistCorrEtaPhi->Fill(deltaPhi, deltaEta, weight);
0ba09863 1921 if (fHistCorrAvgEtaPhi)
1922 fHistCorrAvgEtaPhi->Fill(deltaPhi, avgEta, weight);
ebdf1dbb 1923}
1924
1925void AliAnalysisTaskJetProtonCorr::AliHistCorr::Fill(TLorentzVector *trgPart, AliVParticle *assPart, Float_t weight)
1926{
2751fb9a 1927 // fill correlation histograms from trigger direction and associate particle
1928
ebdf1dbb 1929 Float_t deltaEta = assPart->Eta() - trgPart->Eta();
1930 Float_t avgEta = (assPart->Eta() + trgPart->Eta()) / 2.;
1931 Float_t deltaPhi = assPart->Phi() - trgPart->Phi();
09fd5f9f 1932 while (deltaPhi > (2.*TMath::Pi() + fHistDphiLo))
ebdf1dbb 1933 deltaPhi -= 2. * TMath::Pi();
09fd5f9f 1934 while (deltaPhi < fHistDphiLo)
ebdf1dbb 1935 deltaPhi += 2. * TMath::Pi();
ebdf1dbb 1936
2751fb9a 1937 fHistCorrPhi->Fill(deltaPhi, weight);
0ba09863 1938 if (fHistCorrPhi2) {
1939 Float_t trgPhi = trgPart->Phi();
1940 if (trgPhi < 0)
1941 trgPhi += 2. * TMath::Pi();
1942 fHistCorrPhi2->Fill(trgPhi, assPart->Phi(), weight);
1943 }
ebdf1dbb 1944 fHistCorrEtaPhi->Fill(deltaPhi, deltaEta, weight);
0ba09863 1945 if (fHistCorrAvgEtaPhi)
1946 fHistCorrAvgEtaPhi->Fill(deltaPhi, avgEta, weight);
ebdf1dbb 1947}
1948
1949void AliAnalysisTaskJetProtonCorr::AliHistCorr::Fill(TLorentzVector *trgPart, TLorentzVector *assPart, Float_t weight)
1950{
2751fb9a 1951 // fill correlation histograms from trigger direction and associate direction
1952
ebdf1dbb 1953 Float_t deltaEta = assPart->Eta() - trgPart->Eta();
1954 Float_t avgEta = (assPart->Eta() + trgPart->Eta()) / 2.;
1955 Float_t deltaPhi = assPart->Phi() - trgPart->Phi();
09fd5f9f 1956 if (deltaPhi > (2.*TMath::Pi() + fHistDphiLo))
ebdf1dbb 1957 deltaPhi -= 2. * TMath::Pi();
09fd5f9f 1958 else if (deltaPhi < fHistDphiLo)
ebdf1dbb 1959 deltaPhi += 2. * TMath::Pi();
ebdf1dbb 1960
2751fb9a 1961 fHistCorrPhi->Fill(deltaPhi, weight);
0ba09863 1962 if (fHistCorrPhi2) {
1963 Float_t trgPhi = trgPart->Phi();
1964 if (trgPhi < 0)
1965 trgPhi += 2. * TMath::Pi();
1966 Float_t assPhi = assPart->Phi();
1967 if (assPhi < 0)
1968 assPhi += 2. * TMath::Pi();
1969 fHistCorrPhi2->Fill(trgPhi, assPhi, weight);
1970 }
ebdf1dbb 1971 fHistCorrEtaPhi->Fill(deltaPhi, deltaEta, weight);
0ba09863 1972 if (fHistCorrAvgEtaPhi)
1973 fHistCorrAvgEtaPhi->Fill(deltaPhi, avgEta, weight);
23536a91 1974}
1975
1976// ----- histogram management -----
1977TH1* AliAnalysisTaskJetProtonCorr::AddHistogram(Hist_t hist, const char *hid, TString title,
1978 Int_t xbins, Float_t xmin, Float_t xmax,
1979 Int_t binType)
1980{
1981 TString hName;
1982 hName.Form("%s_%s", fShortTaskId, hid);
1983 hName.ToLower();
1984 TH1 *h = 0x0;
1985 if (binType == 0)
1986 h = new TH1I(hName.Data(), title,
1987 xbins, xmin, xmax);
1988 else
1989 h = new TH1F(hName.Data(), title,
1990 xbins, xmin, xmax);
1991 GetHistogram(hist) = h;
1992 fOutputList->Add(h);
1993 return h;
1994}
1995
1996TH2* AliAnalysisTaskJetProtonCorr::AddHistogram(Hist_t hist, const char *hid, TString title,
1997 Int_t xbins, Float_t xmin, Float_t xmax,
1998 Int_t ybins, Float_t ymin, Float_t ymax,
1999 Int_t binType)
2000{
2001 TString hName;
2002 hName.Form("%s_%s", fShortTaskId, hid);
2003 hName.ToLower();
2004 TH1 *h = 0x0;
2005 if (binType == 0)
2006 h = GetHistogram(hist) = new TH2I(hName.Data(), title,
2007 xbins, xmin, xmax,
2008 ybins, ymin, ymax);
2009 else
2010 h = GetHistogram(hist) = new TH2F(hName.Data(), title,
2011 xbins, xmin, xmax,
2012 ybins, ymin, ymax);
2013 fOutputList->Add(h);
2014 return (TH2*) h;
2015}
2016
2017TH3* AliAnalysisTaskJetProtonCorr::AddHistogram(Hist_t hist, const char *hid, TString title,
2018 Int_t xbins, Float_t xmin, Float_t xmax,
2019 Int_t ybins, Float_t ymin, Float_t ymax,
2020 Int_t zbins, Float_t zmin, Float_t zmax,
2021 Int_t binType)
2022{
2023 TString hName;
2024 hName.Form("%s_%s", fShortTaskId, hid);
2025 hName.ToLower();
2026 TH1 *h = 0x0;
2027 if (binType == 0)
2028 h = GetHistogram(hist) = new TH3I(hName.Data(), title,
2029 xbins, xmin, xmax,
2030 ybins, ymin, ymax,
2031 zbins, zmin, zmax);
2032 else
2033 h = GetHistogram(hist) = new TH3F(hName.Data(), title,
2034 xbins, xmin, xmax,
2035 ybins, ymin, ymax,
2036 zbins, zmin, zmax);
2037 fOutputList->Add(h);
2038 return (TH3*) h;
2039}