Include L0 trigger in the trigger maker and the trigger patch object Add monitoring...
[u/mrichter/AliRoot.git] / PWGJE / UserTasks / AliAnalysisTaskJetProtonCorr.cxx
1 // ROOT
2 #include "TFile.h"
3 #include "TList.h"
4 #include "TH1.h"
5 #include "TH2.h"
6 #include "TH3.h"
7 #include "TF1.h"
8 #include "TFormula.h"
9 #include "TRandom.h"
10 #include "TSpline.h"
11
12 // analysis framework
13 #include "AliAnalysisManager.h"
14 #include "AliInputEventHandler.h"
15 #include "AliVEvent.h"
16 #include "AliVTrack.h"
17 #include "AliVTrdTrack.h"
18 #include "AliVVertex.h"
19 #include "AliPIDResponse.h"
20 #include "AliEventPoolManager.h"
21 #include "AliOADBContainer.h"
22 #include "AliTOFPIDParams.h"
23 #include "AliAnalysisTaskVnV0.h"
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
49 #include <iostream>
50 #include <cmath>
51
52 AliAnalysisTaskJetProtonCorr::AliAnalysisTaskJetProtonCorr(const char *name) :
53   AliAnalysisTaskSE(name),
54   fMCEvent(0x0),
55   fESDEvent(0x0),
56   fAODEvent(0x0),
57   fRunNumber(-1),
58   fOADBContainerTOF(0x0),
59   fParamsTOF(0x0),
60   fEventplane(0x0),
61   fTriggerMask(0),
62   fClassMask(0),
63   fCentrality(100.),
64   fCentralityCheck(100.),
65   fZvtx(0.),
66   fPIDResponse(0x0),
67   fEventPlaneAngle(5.),
68   fEventPlaneRes(0.),
69   fEventPlaneAngleCheck(5.),
70   fEventPlaneAngle3(5.),
71   fPrimTrackArrayAss(0x0),
72   fPrimTrackArrayTrg(0x0),
73   fPrimConstrainedTrackArray(new TClonesArray("AliESDtrack", 100)),
74   fJetArray(0x0),
75   fPoolMgr(),
76   fPool(),
77   fHistCorr(0x0),
78   fErrorMsg(10),
79   fOutputList(),
80   fHist(),
81   fShortTaskId("jet_prot_corr"),
82   fUseStandardCuts(kTRUE),
83   fUseEvplaneV0(kFALSE),
84   fCutsPrimTrg(0x0),
85   fCutsPrimTrgConstrain(new AliESDtrackCuts()),
86   fCutsPrimAss(0x0),
87   fCutsTwoTrackEff(0.02),
88   fAssFilterMask(1 << 10),
89   fRequirePID(kTRUE),
90   fTrgJetEtaMax(0.45),
91   fHadEtaMax(0.8),
92   fTrgPartPtMin(6.),
93   fTrgPartPtMax(8.),
94   fTrgJetPtMin(50.),
95   fTrgJetPtMax(80.),
96   fTrgJetLeadTrkPtMin(6.),
97   fTrgJetLeadTrkPtMax(100.),
98   fTrgJetAreaMin(0.6 * TMath::Pi() * 0.2*0.2),
99   fAssPartPtMin(2.),
100   fAssPartPtMax(4.),
101   fTrgAngleToEvPlane(TMath::Pi() / 4.),
102   fToyMeanNoPart(.5),
103   fToyRadius(.8),
104   fToySmearPhi(.2),
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())),
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)
114 {
115   // default ctor
116
117   fkCorrTypeName[kCorrHadHad]  = "hh";
118   fkCorrTypeName[kCorrHadProt] = "hp";
119   fkCorrTypeName[kCorrJetHad]  = "jh";
120   fkCorrTypeName[kCorrJetProt] = "jp";
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";
129
130   fkClassName[kClCentral]      = "cent";
131   fkClassName[kClSemiCentral]  = "semi";
132   // fkClassName[kClDijet]        = "dijet";
133
134   fkEvName[kEvSame] = "same";
135   fkEvName[kEvMix]  = "mixed";
136
137   // track cuts for associates
138   if (fUseStandardCuts) {
139     fCutsPrimAss = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
140   } else {
141     fCutsPrimAss = new AliESDtrackCuts();
142
143     // this is taken from PWGJE track cuts
144     TFormula *f1NClustersTPCLinearPtDep = new TFormula("f1NClustersTPCLinearPtDep","70.+30./20.*x");
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);
152     // ITS
153     fCutsPrimAss->SetRequireITSRefit(kTRUE);
154     //accept secondaries
155     fCutsPrimAss->SetMaxDCAToVertexXY(2.4);
156     fCutsPrimAss->SetMaxDCAToVertexZ(3.2);
157     fCutsPrimAss->SetDCAToVertex2D(kTRUE);
158     //reject fakes
159     fCutsPrimAss->SetMaxChi2PerClusterITS(36);
160     fCutsPrimAss->SetMaxChi2TPCConstrainedGlobal(36);
161
162     fCutsPrimAss->SetRequireSigmaToVertex(kFALSE);
163
164     fCutsPrimAss->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
165   }
166
167   fCutsPrimAss->SetEtaRange(-fHadEtaMax, fHadEtaMax);
168   fCutsPrimAss->SetPtRange(0.15, 1E+15);
169
170   // track cuts for triggers
171   fCutsPrimTrg = new AliESDtrackCuts(*fCutsPrimAss);
172
173   // azimuthal modulation for triggers
174   fTrgJetPhiModCent->SetParameter(0, fTrgJetV2Cent);
175   fTrgJetPhiModSemi->SetParameter(0, fTrgJetV2Semi);
176   fTrgHadPhiModCent->SetParameter(0, fTrgHadV2Cent);
177   fTrgHadPhiModSemi->SetParameter(0, fTrgHadV2Semi);
178
179   // event mixing pool
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   };
185   Int_t nCentralityBins = sizeof(centralityBins)/sizeof(centralityBins[0]);
186
187   Double_t vertexBins[] = {
188     -10., -8., -6., -4., -2., 0., 2., 4., 6., 8., 10.
189   };
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
197   const Int_t poolSize = 10; // unused by current event pool implementation
198   const Int_t trackDepth = 50000; // number of tracks to maintain in mixing buffer
199   const Float_t trackDepthFraction = 0.1;
200   const Int_t targetEvents = 1;
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   }
208
209   fHistCorr = new AliHistCorr*[kEvLast*kCorrLast*kClLast];
210
211   DefineOutput(1, TList::Class());
212 }
213
214 AliAnalysisTaskJetProtonCorr::~AliAnalysisTaskJetProtonCorr()
215 {
216   // dtor
217
218   // delete [] fHistCorr;
219 }
220
221 void AliAnalysisTaskJetProtonCorr::UserCreateOutputObjects()
222 {
223   // create user output objects
224
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
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));
250   histStat->GetXaxis()->SetBinLabel(ID(kStatVtx));
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));
256   histStat->GetXaxis()->SetBinLabel(ID(kStatCentral));
257   histStat->GetXaxis()->SetBinLabel(ID(kStatSemiCentral));
258
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
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));
271   // hist->GetYaxis()->SetBinLabel(LAB(kClDijet));
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));
279   // hist->GetYaxis()->SetBinLabel(LAB(kClDijet));
280   AddHistogram(ID(kHistCentralityVsMult), "centrality - multiplicity;centrality percentile (%);N_{prim}",
281                100, 0., 100.,
282                100, 0., 2500.);
283
284   AddHistogram(ID(kHistSignalTPC), "TPC dE/dx;p (GeV/#it{c});dE/dx (arb. units)",
285                100, 0., 10., 200, 0., 300.);
286   AddHistogram(ID(kHistSignalTOF), "TOF time of flight;p_{T} (GeV/#it{c});t (ns)",
287                100, 0., 10., 200, 0., 50.);
288   AddHistogram(ID(kHistBetaTOF), "TOF beta;p (GeV/#it{c}); #beta",
289                100, 0., 10.,
290                100, 0., 1.);
291
292   AddHistogram(ID(kHistDeltaTPC), "TPC dE/dx;p (GeV/#it{c});dE/dx (arb. units)",
293                100, 0., 10., 200, -100., 100.);
294   AddHistogram(ID(kHistDeltaTPCSemi), "TPC dE/dx;p (GeV/#it{c});dE/dx (arb. units)",
295                100, 0., 10., 200, -100., 100.);
296   AddHistogram(ID(kHistDeltaTOF), "TOF time of flight;p (GeV/#it{c});t (ns)",
297                100, 0., 10., 200, -2., 2.);
298   AddHistogram(ID(kHistDeltaTOFSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
299                100, 0., 10., 200, -2., 2.);
300
301   // Nsigma templates - central
302   AddHistogram(ID(kHistNsigmaTPCe), "TPC N_{#sigma,p} - e hypothesis;p (GeV/#it{c});N_{#sigma,p}",
303                100, 0., 10.,
304                100, -25., 25.);
305   AddHistogram(ID(kHistNsigmaTPCmu), "TPC N_{#sigma,p} - #mu hypothesis;p (GeV/#it{c});N_{#sigma,p}",
306                100, 0., 10.,
307                100, -25., 25.);
308   AddHistogram(ID(kHistNsigmaTPCpi), "TPC N_{#sigma,p} - #pi hypothesis;p (GeV/#it{c});N_{#sigma,p}",
309                100, 0., 10.,
310                100, -25., 25.);
311   AddHistogram(ID(kHistNsigmaTPCk), "TPC N_{#sigma,p} - K hypothesis;p (GeV/#it{c});N_{#sigma,p}",
312                100, 0., 10.,
313                100, -25., 25.);
314   AddHistogram(ID(kHistNsigmaTPCp), "TPC N_{#sigma,p} - p hypothesis;p (GeV/#it{c});N_{#sigma,p}",
315                100, 0., 10.,
316                100, -25., 25.);
317   AddHistogram(ID(kHistNsigmaTPCd), "TPC N_{#sigma,p} - d hypothesis;p (GeV/#it{c});N_{#sigma,p}",
318                100, 0., 10.,
319                100, -25., 25.);
320   AddHistogram(ID(kHistNsigmaTPCe_e), "TPC N_{#sigma,p} - e hypothesis (id. e);p (GeV/#it{c});N_{#sigma,p}",
321                100, 0., 10.,
322                100, -25., 25.);
323
324   AddHistogram(ID(kHistNsigmaTOFe), "TOF N_{#sigma,p} - e hypothesis;p (GeV/#it{c});N_{#sigma,p}",
325                100, 0., 10.,
326                200, -100., 100.);
327   AddHistogram(ID(kHistNsigmaTOFmu), "TOF N_{#sigma,p} - #mu hypothesis;p (GeV/#it{c});N_{#sigma,p}",
328                100, 0., 10.,
329                200, -100., 100.);
330   AddHistogram(ID(kHistNsigmaTOFpi), "TOF N_{#sigma,p} - #pi hypothesis;p (GeV/#it{c});N_{#sigma,p}",
331                100, 0., 10.,
332                200, -100., 100.);
333   AddHistogram(ID(kHistNsigmaTOFk), "TOF N_{#sigma,p} - K hypothesis;p (GeV/#it{c});N_{#sigma,p}",
334                100, 0., 10.,
335                200, -100., 100.);
336   AddHistogram(ID(kHistNsigmaTOFp), "TOF N_{#sigma,p} - p hypothesis;p (GeV/#it{c});N_{#sigma,p}",
337                100, 0., 10.,
338                200, -100., 100.);
339   AddHistogram(ID(kHistNsigmaTOFd), "TOF N_{#sigma,p} - d hypothesis;p (GeV/#it{c});N_{#sigma,p}",
340                100, 0., 10.,
341                200, -100., 100.);
342   AddHistogram(ID(kHistNsigmaTOFmismatch), "TOF N_{#sigma,p} - mismatch;p (GeV/#it{c});N_{#sigma,p}",
343                100, 0., 10.,
344                200, -100., 100.);
345
346   // Nsigma templates - semi-central
347   AddHistogram(ID(kHistNsigmaTPCeSemi), "TPC N_{#sigma,p} - e hypothesis;p (GeV/#it{c});N_{#sigma,p}",
348                100, 0., 10.,
349                100, -25., 25.);
350   AddHistogram(ID(kHistNsigmaTPCmuSemi), "TPC N_{#sigma,p} - #mu hypothesis;p (GeV/#it{c});N_{#sigma,p}",
351                100, 0., 10.,
352                100, -25., 25.);
353   AddHistogram(ID(kHistNsigmaTPCpiSemi), "TPC N_{#sigma,p} - #pi hypothesis;p (GeV/#it{c});N_{#sigma,p}",
354                100, 0., 10.,
355                100, -25., 25.);
356   AddHistogram(ID(kHistNsigmaTPCkSemi), "TPC N_{#sigma,p} - K hypothesis;p (GeV/#it{c});N_{#sigma,p}",
357                100, 0., 10.,
358                100, -25., 25.);
359   AddHistogram(ID(kHistNsigmaTPCpSemi), "TPC N_{#sigma,p} - p hypothesis;p (GeV/#it{c});N_{#sigma,p}",
360                100, 0., 10.,
361                100, -25., 25.);
362   AddHistogram(ID(kHistNsigmaTPCdSemi), "TPC N_{#sigma,p} - d hypothesis;p (GeV/#it{c});N_{#sigma,p}",
363                100, 0., 10.,
364                100, -25., 25.);
365   AddHistogram(ID(kHistNsigmaTPCe_eSemi), "TPC N_{#sigma,p} - e hypothesis (id. e);p (GeV/#it{c});N_{#sigma,p}",
366                100, 0., 10.,
367                100, -25., 25.);
368
369   AddHistogram(ID(kHistNsigmaTOFeSemi), "TOF N_{#sigma,p} - e hypothesis;p (GeV/#it{c});N_{#sigma,p}",
370                100, 0., 10.,
371                200, -100., 100.);
372   AddHistogram(ID(kHistNsigmaTOFmuSemi), "TOF N_{#sigma,p} - #mu hypothesis;p (GeV/#it{c});N_{#sigma,p}",
373                100, 0., 10.,
374                200, -100., 100.);
375   AddHistogram(ID(kHistNsigmaTOFpiSemi), "TOF N_{#sigma,p} - #pi hypothesis;p (GeV/#it{c});N_{#sigma,p}",
376                100, 0., 10.,
377                200, -100., 100.);
378   AddHistogram(ID(kHistNsigmaTOFkSemi), "TOF N_{#sigma,p} - K hypothesis;p (GeV/#it{c});N_{#sigma,p}",
379                100, 0., 10.,
380                200, -100., 100.);
381   AddHistogram(ID(kHistNsigmaTOFpSemi), "TOF N_{#sigma,p} - p hypothesis;p (GeV/#it{c});N_{#sigma,p}",
382                100, 0., 10.,
383                200, -100., 100.);
384   AddHistogram(ID(kHistNsigmaTOFdSemi), "TOF N_{#sigma,p} - d hypothesis;p (GeV/#it{c});N_{#sigma,p}",
385                100, 0., 10.,
386                200, -100., 100.);
387   AddHistogram(ID(kHistNsigmaTOFmismatchSemi), "TOF N_{#sigma,p} - mismatch;p (GeV/#it{c});N_{#sigma,p}",
388                100, 0., 10.,
389                200, -100., 100.);
390
391   // delta templates
392   AddHistogram(ID(kHistDeltaTOFe), "TOF #Delta;p (GeV/#it{c});t (ns)",
393                100, 0., 10., 200, -2., 2.);
394   AddHistogram(ID(kHistDeltaTOFmu), "TOF #Delta;p (GeV/#it{c});t (ns)",
395                100, 0., 10., 200, -2., 2.);
396   AddHistogram(ID(kHistDeltaTOFpi), "TOF #Delta;p (GeV/#it{c});t (ns)",
397                100, 0., 10., 200, -2., 2.);
398   AddHistogram(ID(kHistDeltaTOFk), "TOF #Delta;p (GeV/#it{c});t (ns)",
399                100, 0., 10., 200, -2., 2.);
400   AddHistogram(ID(kHistDeltaTOFp), "TOF #Delta;p (GeV/#it{c});t (ns)",
401                100, 0., 10., 200, -2., 2.);
402   AddHistogram(ID(kHistDeltaTOFd), "TOF #Delta;p (GeV/#it{c});t (ns)",
403                100, 0., 10., 200, -2., 2.);
404
405   AddHistogram(ID(kHistDeltaTOFeSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
406                100, 0., 10., 200, -2., 2.);
407   AddHistogram(ID(kHistDeltaTOFmuSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
408                100, 0., 10., 200, -2., 2.);
409   AddHistogram(ID(kHistDeltaTOFpiSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
410                100, 0., 10., 200, -2., 2.);
411   AddHistogram(ID(kHistDeltaTOFkSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
412                100, 0., 10., 200, -2., 2.);
413   AddHistogram(ID(kHistDeltaTOFpSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
414                100, 0., 10., 200, -2., 2.);
415   AddHistogram(ID(kHistDeltaTOFdSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
416                100, 0., 10., 200, -2., 2.);
417
418   // sigma comparisons
419   // AddHistogram(ID(kHistExpSigmaTOFe), "TOF time of flight;p (GeV/#it{c});t (ns)",
420   //           100, 0., 10., 200, 0., .25);
421   // AddHistogram(ID(kHistExpSigmaTOFmu), "TOF time of flight;p (GeV/#it{c});t (ns)",
422   //           100, 0., 10., 200, 0., .25);
423   // AddHistogram(ID(kHistExpSigmaTOFpi), "TOF time of flight;p (GeV/#it{c});t (ns)",
424   //           100, 0., 10., 200, 0., .25);
425   // AddHistogram(ID(kHistExpSigmaTOFk), "TOF time of flight;p (GeV/#it{c});t (ns)",
426   //           100, 0., 10., 200, 0., .25);
427   // AddHistogram(ID(kHistExpSigmaTOFp), "TOF time of flight;p (GeV/#it{c});t (ns)",
428   //           100, 0., 10., 200, 0., .25);
429   // AddHistogram(ID(kHistExpSigmaTOFd), "TOF time of flight;p (GeV/#it{c});t (ns)",
430   //           100, 0., 10., 200, 0., .25);
431
432   // AddHistogram(ID(kHistExpSigmaTOFeSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
433   //           100, 0., 10., 200, 0., .25);
434   // AddHistogram(ID(kHistExpSigmaTOFmuSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
435   //           100, 0., 10., 200, 0., .25);
436   // AddHistogram(ID(kHistExpSigmaTOFpiSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
437   //           100, 0., 10., 200, 0., .25);
438   // AddHistogram(ID(kHistExpSigmaTOFkSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
439   //           100, 0., 10., 200, 0., .25);
440   // AddHistogram(ID(kHistExpSigmaTOFpSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
441   //           100, 0., 10., 200, 0., .25);
442   // AddHistogram(ID(kHistExpSigmaTOFdSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
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);
470
471   // Nsigma distributions
472   AddHistogram(ID(kHistNsigmaTPCTOF), "N_{#sigma,p} TPC-TOF;p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
473                100, 0., 10.,
474                100, -25., 25.,
475                200, -100., 100.);
476   AddHistogram(ID(kHistNsigmaTPCTOFPt), "N_{#sigma,p} TPC-TOF;p_{T} (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
477                100, 0., 10.,
478                100, -25., 25.,
479                200, -100., 100.);
480   AddHistogram(ID(kHistNsigmaTPCTOFUsed), "N_{#sigma,p} TPC-TOF;p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
481                100, 0., 10.,
482                100, -25., 25.,
483                200, -100., 100.);
484   AddHistogram(ID(kHistNsigmaTPCTOFUsedCentral), "N_{#sigma,p} TPC-TOF (central);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
485                100, 0., 10.,
486                100, -25., 25.,
487                200, -100., 100.);
488   AddHistogram(ID(kHistNsigmaTPCTOFUsedSemiCentral), "N_{#sigma,p} TPC-TOF (semi-central);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
489                100, 0., 10.,
490                100, -25., 25.,
491                200, -100., 100.);
492   AddHistogram(ID(kHistNsigmaTPCTOFUsedPt), "N_{#sigma,p} TPC-TOF;p_{T} (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
493                50, 0., 10.,
494                100, -25., 25.,
495                200, -100., 100.);
496   AddHistogram(ID(kHistNsigmaTPCTOFUsedPtCentral), "N_{#sigma,p} TPC-TOF;p_{T} (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
497                50, 0., 10.,
498                100, -25., 25.,
499                200, -100., 100.);
500   AddHistogram(ID(kHistNsigmaTPCTOFUsedPtSemiCentral), "N_{#sigma,p} TPC-TOF;p_{T} (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
501                50, 0., 10.,
502                100, -25., 25.,
503                200, -100., 100.);
504
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
555   AddHistogram(ID(kHistNevMix), "no. of events for mixing;N_{ev};counts",
556                20, -.5, 19.5,
557                2, -.5, 1.5);
558
559   AddHistogram(ID(kHistEvPlane), "default event plane;#Psi;counts",
560                100, -0. * TMath::Pi(), 1. * TMath::Pi());
561   AddHistogram(ID(kHistEvPlaneRes), "resolution of default event plane;R^{2};centrality (%)",
562                100, -1., 1.,
563                50, 0., 100.);
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);
572   AddHistogram(ID(kHistEvPlane3), "3rd order event plane;#Psi;counts",
573                100, -0. * TMath::Pi(), 2./3. * TMath::Pi());
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(),
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.);
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",
586                100, -0. * TMath::Pi(), 1. * TMath::Pi(),
587                100, -0. * TMath::Pi(), 1. * TMath::Pi(),
588                kClLast, -.5, kClLast-.5);
589
590   AddHistogram(ID(kHistJetPtCentral), "jet spectrum - central;p_{T}^{jet,ch} (GeV/#it{c});counts",
591                40, 0., 200.);
592   AddHistogram(ID(kHistJetPtSemi), "jet spectrum - semi-peripheral;p_{T}^{jet,ch} (GeV/#it{c});counts",
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
608   AddHistogram(ID(kHistPhiTrgJetEvPlane), "trg jet;#varphi - #Psi_{ev};centrality",
609                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
610                100, 0., 100.);
611   AddHistogram(ID(kHistPhiTrgHadEvPlane), "trg had;#varphi - #Psi_{ev};centrality",
612                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
613                100, 0., 100.);
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.);
620   AddHistogram(ID(kHistPhiAssHadEvPlane), "ass had;#varphi - #Psi_{ev};centrality",
621                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
622                100, 0., 100.);
623   AddHistogram(ID(kHistPhiAssProtEvPlane), "ass prot;#varphi - #Psi_{ev};centrality",
624                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
625                100, 0., 100.);
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(),
629                100, 0., 100.);
630
631   AddHistogram(ID(kHistPhiTrgJetEvPlane3), "trg jet;#varphi - #Psi_{ev};centrality",
632                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
633                100, 0., 100.);
634   AddHistogram(ID(kHistPhiTrgHadEvPlane3), "trg had;#varphi - #Psi_{ev};centrality",
635                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
636                100, 0., 100.);
637   AddHistogram(ID(kHistPhiAssHadEvPlane3), "ass had;#varphi - #Psi_{ev};centrality",
638                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
639                100, 0., 100.);
640   AddHistogram(ID(kHistPhiAssProtEvPlane3), "ass prot;#varphi - #Psi_{ev};centrality",
641                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
642                100, 0., 100.);
643
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) {
647         // we don't need the mixed event histograms for the embedded excess particles
648         if ((iCorr > kCorrRndJetProt) && (iEv == kEvMix))
649           continue;
650
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
660 Bool_t AliAnalysisTaskJetProtonCorr::Notify()
661 {
662   // actions to be taken upon notification about input file change
663
664   return AliAnalysisTaskSE::Notify();
665 }
666
667 void 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
695 void 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
734     FillH1(kHistVertexZ, fZvtx);
735
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();
746     if (IsClass(kClCentral))
747       FillH1(kHistStat, kStatCentral);
748     if (IsClass(kClSemiCentral))
749       FillH1(kHistStat, kStatSemiCentral);
750
751     FillH1(kHistEvPlane, fEventPlaneAngle);
752     FillH2(kHistEvPlaneRes, fEventPlaneRes, fCentrality);
753     FillH1(kHistEvPlaneCheck, fEventPlaneAngleCheck);
754     FillH1(kHistEvPlane3, fEventPlaneAngle3);
755     FillH2(kHistEvPlaneCross, TMath::Cos(2. * (fEventPlaneAngle - fEventPlaneAngleCheck)), GetCentrality());
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);
760         FillH2(kHistEvPlaneUsed, fEventPlaneAngle, iClass);
761         FillH2(kHistEvPlaneCheckUsed, fEventPlaneAngleCheck, iClass);
762         FillH3(kHistEvPlaneCorr, fEventPlaneAngle, fEventPlaneAngleCheck, iClass);
763       }
764     }
765
766     Bool_t corrEta  = fPIDResponse->UseTPCEtaCorrection();
767     Bool_t corrMult = fPIDResponse->UseTPCMultiplicityCorrection();
768
769     // select trigger particles and potential associated particles/protons
770     TObjArray trgArray[kTrgLast];
771     TObjArray assArray[kAssLast];
772
773     // prepare TOF response
774     TF1 fTOFsignal("fTOFsignal", &AliAnalysisTaskJetProtonCorr::TOFsignal, -2440., 2440., 4);
775     fTOFsignal.SetParameter(0, 1.);
776     fTOFsignal.SetParameter(1, 0.);
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     }
787
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);
793       FillH2(kHistSignalTPC, trk->P(), trk->GetTPCsignal());
794       // ??? pt or p?
795       FillH2(kHistSignalTOF, trk->P(), trk->GetTOFsignal() * 1.e-3); // ps -> ns
796       FillH3(kHistNsigmaTPCTOF,
797              trk->P(),
798              fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
799              fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
800       FillH3(kHistNsigmaTPCTOFPt,
801              trk->Pt(),
802              fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
803              fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
804
805       FillH3(kHistPhiAssHadVsEvPlane, fEventPlaneAngle, trk->Phi(), fCentrality);
806
807       if (AcceptAssoc(trk)) {
808         Float_t phiRel = GetPhiRel2(trk);
809         Float_t phiRel3 = trk->Phi() - fEventPlaneAngle3;
810         if (phiRel3 < 0.)
811           phiRel3 += 2. * TMath::Pi();
812         assArray[kAssHad].Add(trk);
813         FillH1(kHistEtaPhiAssHad, trk->Phi(), trk->Eta());
814         FillH2(kHistPhiAssHadEvPlane, phiRel, fCentrality);
815         FillH2(kHistPhiAssHadEvPlane3, phiRel3, fCentrality);
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));
824
825         // central events
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));
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           }
847         }
848
849         // semi-central events
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));
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           }
871         }
872
873         // protons
874         if (IsProton(trk)) {
875           assArray[kAssProt].Add(trk);
876           FillH1(kHistEtaPhiAssProt, trk->Phi(), trk->Eta());
877           FillH2(kHistPhiAssProtEvPlane, phiRel, fCentrality);
878           FillH2(kHistPhiAssProtEvPlane3, phiRel3, fCentrality);
879         }
880
881         // template generation
882         if (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC, trk) == AliPIDResponse::kDetPidOk) {
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
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);
898
899             if(IsClass(kClCentral))
900               FillH2(kHistNsigmaTPCe, trk->P(), (rndTPCx-expTPC)/expSigmaTPC, 1., iParticle);
901             if(IsClass(kClSemiCentral))
902               FillH2(kHistNsigmaTPCeSemi, trk->P(), (rndTPCx-expTPC)/expSigmaTPC, 1., iParticle);
903           }
904         }
905
906         if (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, trk) == AliPIDResponse::kDetPidOk) {
907           Double_t deltaTOF;
908           fPIDResponse->GetSignalDelta(AliPIDResponse::kTOF, trk, AliPID::kProton, deltaTOF);
909
910           AliTOFPIDResponse &tofResponse = fPIDResponse->GetTOFResponse();
911
912           Float_t p = trk->P();
913           if (const AliExternalTrackParam *param = trk->GetInnerParam())
914             p = param->GetP();
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
925           // Double_t signalSigma = TMath::Sqrt(fTOFsignal.Variance(-2000., 2000.));
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);
930           // tracking smearing (default parameters, should be overwritten from OADB)
931           Double_t fPar[] = { 0.008, 0.008, 0.002, 40.0 };
932           if (fParamsTOF)
933             for (Int_t i = 0; i < 4; ++i)
934               fPar[i] = fParamsTOF->GetSigParams(i);
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));
939             // Double_t cmpSigmaTOFx = fPIDResponse->GetTOFResponse().GetExpectedSigma(p, expTOFx, AliPID::EParticleType(iParticle));
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));
946             // Double_t tmpSigmaTOFx = TMath::Sqrt(signalSigma*signalSigma +
947             //                                  timezeroSigma*timezeroSigma +
948             //                                  expSigmaTOFx*expSigmaTOFx + fPar[3]*fPar[3]/p/p);
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);
963               // FillH2(kHistExpSigmaTOFe, p, cmpSigmaTOFx * 1.e-3, 1., iParticle);
964               // FillH2(kHistCmpSigmaTOFe, cmpSigmaTOFx * 1.e-3, tmpSigmaTOFx * 1.e-3, 1., iParticle);
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);
969               // FillH2(kHistCmpSigmaTOFeSemi, cmpSigmaTOFx * 1.e-3, tmpSigmaTOFx * 1.e-3, 1., iParticle);
970               // FillH2(kHistExpSigmaTOFeSemi, p, cmpSigmaTOFx * 1.e-3, 1., iParticle);
971             }
972           }
973
974           Double_t rndTOFmismatch = AliTOFPIDResponse::GetMismatchRandomValue(trk->Eta());
975
976           if(IsClass(kClCentral)) {
977             FillH2(kHistNsigmaTOFmismatch, trk->P(), (rndTOFmismatch - expTOF) / expSigmaTOF);
978             FillH2(kHistDeltaTOF, trk->P(), deltaTOF * 1.e-3); // ps -> ns
979           }
980           if(IsClass(kClSemiCentral)) {
981             FillH2(kHistNsigmaTOFmismatchSemi, trk->P(), (rndTOFmismatch - expTOF) / expSigmaTOF);
982             FillH2(kHistDeltaTOFSemi, trk->P(), deltaTOF * 1.e-3); // ps -> ns
983           }
984         }
985       }
986     }
987
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
997     // select trigger jet
998     // and remove them from the Q vector
999     Int_t nJets = fJetArray ? fJetArray->GetEntries() : 0;
1000     const TVector2 *qVectorOrig = fEventplane->GetQVector();
1001     TVector2 qVector;
1002     if (qVectorOrig)
1003       qVector = *qVectorOrig;
1004
1005     for (Int_t iJet = 0; iJet < nJets; ++iJet) {
1006       AliAODJet *jet = (AliAODJet*) fJetArray->At(iJet);
1007
1008       if (AcceptTrigger(jet)) {
1009         trgArray[kTrgJet].Add(jet);
1010         FillH1(kHistEtaPhiTrgJet, jet->Phi(), jet->Eta());
1011
1012         if (qVectorOrig) {
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             
1017             if (fEventplane && track &&
1018                 (track->GetID() > -1)) {
1019               TVector2 evplaneContrib(fEventplane->GetQContributionX(track),
1020                                       fEventplane->GetQContributionY(track));
1021               qVector -= evplaneContrib;
1022             }
1023           }
1024         }
1025       }
1026     }
1027     if (qVectorOrig) {
1028       for (Int_t iClass = 0; iClass < kClLast; ++iClass) {
1029         if (IsClass((Class_t) iClass)) {
1030           FillH3(kHistEvPlaneCorrNoTrgJets, qVectorOrig->Phi()/2., qVector.Phi()/2., iClass);
1031           if (trgArray[kTrgJet].GetEntriesFast() > 0)
1032             FillH3(kHistEvPlaneCorrNoTrgJetsTrgd, qVectorOrig->Phi()/2., qVector.Phi()/2., iClass);
1033         }
1034       }
1035     }
1036
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
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) {
1050           for (Int_t iAss = 0; iAss <= kAssProt; ++iAss) {
1051             // same event
1052             Correlate((Trg_t) iTrg, (Ass_t) iAss, (Class_t) iClass, kEvSame, &trgArray[iTrg], &assArray[iAss]);
1053
1054             // mixed event
1055             AliEventPool *pool = GetPool((Ass_t) iAss);
1056             if (pool)
1057               FillH2(kHistNevMix, pool->GetCurrentNEvents(), iAss);
1058             if (pool && pool->IsReady()) {
1059               AliDebug(1, Form("----- using pool: %i %i %i -----", iClass, iTrg, iAss));
1060               const Int_t nEvents = pool->GetCurrentNEvents();
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             }
1066           }
1067         }
1068
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();
1080           }
1081         }
1082         // }
1083
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]);
1089       }
1090     }
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();
1098   }
1099
1100  stop:
1101   CleanUpEvent();
1102
1103   PostData(1, fOutputList);
1104 }
1105
1106 void AliAnalysisTaskJetProtonCorr::Terminate(const Option_t * /* option */)
1107 {
1108   // actions at task termination
1109
1110 }
1111
1112 void 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
1119 Bool_t AliAnalysisTaskJetProtonCorr::DetectTriggers()
1120 {
1121   fTriggerMask = 0;
1122
1123   AliVEvent::EOfflineTriggerTypes physSel =
1124     (AliVEvent::EOfflineTriggerTypes) ((AliInputEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1125   // TString trgClasses = InputEvent()->GetFiredTriggerClasses();
1126
1127   // physics selection
1128   if (physSel & (AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral))
1129     MarkTrigger(kTriggerInt);
1130
1131   return kTRUE;
1132 }
1133
1134 Bool_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
1147 Bool_t AliAnalysisTaskJetProtonCorr::PrepareEvent()
1148 {
1149   Bool_t eventGood = kTRUE;
1150
1151   // check for run change
1152   if (fRunNumber != InputEvent()->GetRunNumber()) {
1153     fRunNumber = InputEvent()->GetRunNumber();
1154     SetParamsTOF();
1155   }
1156
1157   // retrieve z-vertex position
1158   fZvtx = 100.;
1159   const AliVVertex *vtx = InputEvent()->GetPrimaryVertex();
1160   if (vtx) {
1161     FillH1(kHistStat, kStatVtx);
1162     FillH1(kHistVertexNctb, vtx->GetNContributors());
1163     if (vtx->GetNContributors() >= 3.)
1164       fZvtx = vtx->GetZ();
1165   }
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   }
1180   else
1181     eventGood = kFALSE;
1182
1183   // retrieve event plane
1184   fEventplane = InputEvent()->GetEventplane();
1185   if (fEventplane) {
1186     fEventPlaneAngle3 = fEventplane->GetEventplane("V0", InputEvent(), 3);
1187
1188     if (fUseEvplaneV0) {
1189       fEventPlaneAngle = fEventplane->GetEventplane("V0", InputEvent());
1190       fEventPlaneAngleCheck = fEventplane->GetEventplane("Q");
1191     }
1192     else {
1193       fEventPlaneAngle = fEventplane->GetEventplane("Q");
1194       fEventPlaneRes   = TMath::Cos(2 * fEventplane->GetQsubRes());
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();
1208     if (fEventPlaneAngleCheck < 0)
1209       fEventPlaneAngleCheck += TMath::Pi();
1210
1211     FillH1(kHistStat, kStatEvPlane);
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) {
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     }
1242   }
1243   else if (fAODEvent) {
1244     // associate candidates
1245     fPrimTrackArrayAss = new TObjArray();
1246     const Int_t nTracksAODAss = fAODEvent->GetNumberOfTracks();
1247     for (Int_t iTrack = 0; iTrack < nTracksAODAss; ++iTrack) {
1248       AliAODTrack *trk = fAODEvent->GetTrack(iTrack);
1249       if (trk->TestFilterBit(fAssFilterMask))
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);
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) {
1269       // still try the output event
1270       if (AODEvent())
1271         fJetArray = dynamic_cast<TClonesArray*> (AODEvent()->FindListObject(fJetBranchName));
1272       if (!fJetArray) {
1273         if (fErrorMsg > 0) {
1274           AliError(Form("no jet branch \"%s\" found, in the AODs are:", fJetBranchName));
1275           --fErrorMsg;
1276         }
1277         if (fDebug > 0) {
1278           fAODEvent->GetList()->Print();
1279           if (AODEvent())
1280             AODEvent()->GetList()->Print();
1281         }
1282       }
1283     }
1284   }
1285
1286   // retrieve event pool for the event category
1287   if (eventGood) {
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;
1293     }
1294   }
1295
1296   return eventGood;
1297 }
1298
1299 Float_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
1318   // ??? protect against negative ID
1319   // but should be handled properly by event plane
1320   if (track && (track->GetID() > -1)) {
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             
1330       if (contrib && (contrib->GetID() > -1)) {
1331         TVector2 evplaneContrib(fEventplane->GetQContributionX(contrib),
1332                                 fEventplane->GetQContributionY(contrib));
1333         qVector -= evplaneContrib;
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
1350 Bool_t AliAnalysisTaskJetProtonCorr::AcceptTrigger(AliVTrack *trg)
1351 {
1352   if (TMath::Abs(trg->Eta()) > fHadEtaMax)
1353     return kFALSE;
1354
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;
1367     // if (gRandom->Rndm() > .5)
1368     //   phiRel -= TMath::Pi();
1369     if (phiRel < 0.)
1370       phiRel += 2. * TMath::Pi();
1371     Float_t phiRel3 = trg->Phi() - fEventPlaneAngle3;
1372     if (phiRel3 < 0.)
1373       phiRel3 += 2. * TMath::Pi();
1374     FillH2(kHistPhiTrgHadEvPlane, GetPhiRel2(trg), fCentrality);
1375     FillH2(kHistPhiTrgHadEvPlane3, phiRel3, fCentrality);
1376   }
1377
1378   return (acceptPt && acceptOrientation);
1379 }
1380
1381 AliVTrack* AliAnalysisTaskJetProtonCorr::GetLeadingTrack(const AliAODJet *jet) const
1382 {
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;
1395     }
1396   }
1397
1398   // retrieve the leading track
1399   return (AliAODTrack*) jet->GetRefTracks()->At(iLeadingTrack);
1400 }
1401
1402 Bool_t AliAnalysisTaskJetProtonCorr::AcceptTrigger(AliAODJet *trg)
1403 {
1404   // restrict eta
1405   if (TMath::Abs(trg->Eta()) > fTrgJetEtaMax)
1406     return kFALSE;
1407
1408   // reject too small jets
1409   if (trg->EffectiveAreaCharged() < fTrgJetAreaMin)
1410     return kFALSE;
1411
1412   // require hard leading track
1413   // (leading track biased jet sample)
1414   // but reject jets with too high pt constituents
1415   const Float_t ptLeadTrack = GetLeadingTrack(trg)->Pt();
1416   if ((ptLeadTrack < fTrgJetLeadTrkPtMin) ||
1417       (ptLeadTrack > fTrgJetLeadTrkPtMax))
1418     return kFALSE;
1419
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;
1425
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;
1435     // if (gRandom->Rndm() > .5)
1436     //   phiRel -= TMath::Pi();
1437     if (phiRel < 0.)
1438       phiRel += 2. * TMath::Pi();
1439     Float_t phiRel3 = trg->Phi() - fEventPlaneAngle3;
1440     if (phiRel3 < 0.)
1441       phiRel3 += 2. * TMath::Pi();
1442     FillH2(kHistPhiTrgJetEvPlane, GetPhiRel2(trg), fCentrality);
1443     FillH2(kHistPhiTrgJetEvPlane3, phiRel3, fCentrality);
1444   }
1445
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   }
1453
1454   return (acceptPt && acceptOrientation);
1455 }
1456
1457 Bool_t AliAnalysisTaskJetProtonCorr::AcceptAssoc(const AliVTrack *trk) const
1458 {
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
1469       return kTRUE;
1470   else
1471     return kFALSE;
1472 }
1473
1474 Bool_t AliAnalysisTaskJetProtonCorr::IsProton(const AliVTrack *trk) const
1475 {
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
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
1491 Bool_t AliAnalysisTaskJetProtonCorr::AcceptAngleToEvPlane(Float_t phi, Float_t psi) const
1492 {
1493   Float_t deltaPhi = phi - psi;
1494   while (deltaPhi < 0.)
1495     deltaPhi += 2 * TMath::Pi();
1496
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);
1500
1501   if (TMath::Abs(deltaPhi) < fTrgAngleToEvPlane)
1502     return kTRUE;
1503   else
1504     return kFALSE;
1505 }
1506
1507 Float_t AliAnalysisTaskJetProtonCorr::GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1,
1508                                                   Float_t phi2, Float_t pt2, Float_t charge2,
1509                                                   Float_t radius, Float_t bSign) const
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
1532 Bool_t AliAnalysisTaskJetProtonCorr::AcceptTwoTracks(const AliVParticle *trgPart, const AliVParticle *assPart) const
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
1544   Float_t deta = trgPart->Eta() - assPart->Eta();
1545
1546   Float_t bSign = (InputEvent()->GetMagneticField() > 0) ? 1 : -1;
1547
1548   // optimization
1549   if (TMath::Abs(deta) < fCutsTwoTrackEff) {
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) {
1562         Float_t dphistarabs = TMath::Abs(GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign));
1563         if (dphistarabs < dphistarminabs)
1564           dphistarminabs = dphistarabs;
1565       }
1566
1567       if (dphistarminabs < fCutsTwoTrackEff)
1568         return kFALSE;
1569     }
1570   }
1571
1572   return kTRUE;
1573 }
1574
1575 Bool_t AliAnalysisTaskJetProtonCorr::Correlate(CorrType_t corr, Class_t cl, Ev_t ev,
1576                                                TCollection *trgArray, TCollection *assArray, Float_t weight)
1577 {
1578   AliHistCorr *histCorr = GetHistCorr(corr, cl, ev);
1579   if (!histCorr) {
1580     AliError(Form("no correlation histograms for corr %i, cl %i, ev %i",
1581                   corr, cl, ev));
1582     return kFALSE;
1583   }
1584
1585   TIter assIter(assArray);
1586   while (TObject *assObj = assIter()) {
1587     if (dynamic_cast<AliVParticle*> (assObj)) {
1588       AliVParticle *assPart = (AliVParticle*) assObj;
1589       histCorr->Ass(assPart->Phi(), assPart->Eta(),
1590                     assPart->Pt() > 0. ? (assPart->Charge() / 3.) / assPart->Pt() : 1.e-4,
1591                     weight);
1592     }
1593     else if (dynamic_cast<TLorentzVector*> (assObj)) {
1594       TLorentzVector *assVec = (TLorentzVector*) assObj;
1595       histCorr->Ass(TVector2::Phi_0_2pi(assVec->Phi()), assVec->Eta(), 0., weight);
1596     }
1597   }
1598
1599   TIter trgIter(trgArray);
1600   while (TObject *trgObj = trgIter()) {
1601     if (AliVParticle *trgPart = dynamic_cast<AliVParticle*> (trgObj)) {
1602       // count the trigger
1603       histCorr->Trigger(trgPart->Phi(), trgPart->Eta(),
1604                         trgPart->Pt() > 0. ? (trgPart->Charge() / 3.) / trgPart->Pt() : 1.e-4,
1605                         weight);
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
1621       histCorr->Trigger(TVector2::Phi_0_2pi(trgVec->Phi()), trgVec->Eta(),
1622                         0.,
1623                         weight);
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       }
1635     }
1636   }
1637
1638   return kTRUE;
1639 }
1640
1641 Bool_t AliAnalysisTaskJetProtonCorr::Correlate(Trg_t trg, Ass_t ass, Class_t cl, Ev_t ev,
1642                                                TCollection *trgArray, TCollection *assArray, Float_t weight)
1643 {
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
1650   CorrType_t corr = (CorrType_t) (2 * trg + ass);
1651
1652   return Correlate(corr, cl, ev, trgArray, assArray, weight);
1653 }
1654
1655 Bool_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 {
1663   const Int_t nPart = gRandom->Poisson(fToyMeanNoPart);
1664
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
1681   Float_t trgJetEta = gRandom->Uniform(-fTrgJetEtaMax, fTrgJetEtaMax);
1682   Float_t trgJetPhi = 0.;
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()));
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
1707   Float_t trgJetEtaAway = gRandom->Uniform(-fHadEtaMax, fHadEtaMax);
1708   Float_t trgJetPhiAway = trgJetPhi + TMath::Pi() + gRandom->Gaus(0., fToySmearPhi);
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;
1714     const Float_t r = fToyRadius;
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;
1723     if (TMath::Abs(eta) > fHadEtaMax) {
1724       delete assPart;
1725       continue;
1726     }
1727     while (phi < 0.)
1728       phi += 2. * TMath::Pi();
1729     while (phi > 2. * TMath::Pi())
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
1738   Float_t trgHadEta = gRandom->Uniform(-fHadEtaMax, fHadEtaMax);
1739   Float_t trgHadPhi = 0.;
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()));
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
1764   Float_t trgHadEtaAway = gRandom->Uniform(-fHadEtaMax, fHadEtaMax);
1765   Float_t trgHadPhiAway = trgHadPhi + TMath::Pi() + gRandom->Gaus(0., fToySmearPhi);
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;
1771     const Float_t r = fToyRadius;
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;
1780     if (TMath::Abs(eta) > fHadEtaMax) {
1781       delete assPart;
1782       continue;
1783     }
1784     while (phi < 0.)
1785       phi += 2. * TMath::Pi();
1786     while (phi > 2. * TMath::Pi())
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
1797 Bool_t AliAnalysisTaskJetProtonCorr::CleanUpEvent()
1798 {
1799   if (fAODEvent) {
1800     delete fPrimTrackArrayAss;
1801     fPrimTrackArrayAss = 0x0;
1802     delete fPrimTrackArrayTrg;
1803     fPrimTrackArrayTrg = 0x0;
1804   }
1805
1806   fPrimConstrainedTrackArray->Delete();
1807
1808   return kTRUE;
1809 }
1810
1811 TObjArray* 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++) {
1818     // WARNING: TObject::Clone() is very!!! expensive, unusable
1819     // tracksClone->Add(particle->Clone());
1820
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));
1828   }
1829
1830   return tracksClone;
1831 }
1832
1833 AliAnalysisTaskJetProtonCorr::AliHistCorr::AliHistCorr(TString name, TList *outputList) :
1834   TNamed(name, name),
1835   fOutputList(outputList),
1836   fHistStat(0x0),
1837   fHistCorrPhi(0x0),
1838   fHistCorrPhi2(0x0),
1839   fHistCorrEtaPhi(0x0),
1840   fHistCorrAvgEtaPhi(0x0),
1841   fHistCorrTrgEtaPhi(0x0),
1842   fHistCorrAssEtaPhi(0x0),
1843   fHistCorrTrgEtaPhiQpt(0x0),
1844   fHistCorrAssEtaPhiQpt(0x0),
1845   fHistDphiLo(-TMath::Pi() / 2.),
1846   fHistDphiNbins(120),
1847   fHistDetaNbins(120)
1848 {
1849   // ctor
1850
1851   fHistStat = new TH1F(Form("%s_stat", name.Data()), "statistics",
1852                        1, .5, 1.5);
1853
1854   fHistCorrPhi = new TH1F(Form("%s_phi", name.Data()), ";#Delta#varphi",
1855                           fHistDphiNbins, fHistDphiLo, 2.*TMath::Pi() + fHistDphiLo);
1856   fHistCorrPhi->Sumw2();
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",
1862                              fHistDphiNbins, fHistDphiLo, 2*TMath::Pi() + fHistDphiLo,
1863                              fHistDetaNbins, -2., 2.);
1864   fHistCorrEtaPhi->Sumw2();
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();
1888
1889   fOutputList->Add(fHistStat);
1890   fOutputList->Add(fHistCorrPhi);
1891   fOutputList->Add(fHistCorrPhi2);
1892   fOutputList->Add(fHistCorrEtaPhi);
1893   fOutputList->Add(fHistCorrAvgEtaPhi);
1894   fOutputList->Add(fHistCorrTrgEtaPhi);
1895   fOutputList->Add(fHistCorrAssEtaPhi);
1896   // fOutputList->Add(fHistCorrTrgEtaPhiQpt);
1897   // fOutputList->Add(fHistCorrAssEtaPhiQpt);
1898 }
1899
1900 AliAnalysisTaskJetProtonCorr::AliHistCorr::~AliHistCorr()
1901 {
1902   // dtor
1903 }
1904
1905 void AliAnalysisTaskJetProtonCorr::AliHistCorr::Fill(AliVParticle *trgPart, AliVParticle *assPart, Float_t weight)
1906 {
1907   // fill correlation histograms from trigger particle and associate particle
1908
1909   Float_t deltaEta = assPart->Eta() - trgPart->Eta();
1910   Float_t avgEta = (assPart->Eta() + trgPart->Eta()) / 2.;
1911   Float_t deltaPhi = assPart->Phi() - trgPart->Phi();
1912   while (deltaPhi > (2.*TMath::Pi() + fHistDphiLo))
1913     deltaPhi -= 2. * TMath::Pi();
1914   while (deltaPhi < fHistDphiLo)
1915     deltaPhi += 2. * TMath::Pi();
1916
1917   fHistCorrPhi->Fill(deltaPhi, weight);
1918   if (fHistCorrPhi2)
1919     fHistCorrPhi2->Fill(trgPart->Phi(), assPart->Phi(), weight);
1920   fHistCorrEtaPhi->Fill(deltaPhi, deltaEta, weight);
1921   if (fHistCorrAvgEtaPhi)
1922     fHistCorrAvgEtaPhi->Fill(deltaPhi, avgEta, weight);
1923 }
1924
1925 void AliAnalysisTaskJetProtonCorr::AliHistCorr::Fill(TLorentzVector *trgPart, AliVParticle *assPart, Float_t weight)
1926 {
1927   // fill correlation histograms from trigger direction and associate particle
1928
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();
1932   while (deltaPhi > (2.*TMath::Pi() + fHistDphiLo))
1933     deltaPhi -= 2. * TMath::Pi();
1934   while (deltaPhi < fHistDphiLo)
1935     deltaPhi += 2. * TMath::Pi();
1936
1937   fHistCorrPhi->Fill(deltaPhi, weight);
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   }
1944   fHistCorrEtaPhi->Fill(deltaPhi, deltaEta, weight);
1945   if (fHistCorrAvgEtaPhi)
1946     fHistCorrAvgEtaPhi->Fill(deltaPhi, avgEta, weight);
1947 }
1948
1949 void AliAnalysisTaskJetProtonCorr::AliHistCorr::Fill(TLorentzVector *trgPart, TLorentzVector *assPart, Float_t weight)
1950 {
1951   // fill correlation histograms from trigger direction and associate direction
1952
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();
1956   if (deltaPhi > (2.*TMath::Pi() + fHistDphiLo))
1957     deltaPhi -= 2. * TMath::Pi();
1958   else if (deltaPhi < fHistDphiLo)
1959     deltaPhi += 2. * TMath::Pi();
1960
1961   fHistCorrPhi->Fill(deltaPhi, weight);
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   }
1971   fHistCorrEtaPhi->Fill(deltaPhi, deltaEta, weight);
1972   if (fHistCorrAvgEtaPhi)
1973     fHistCorrAvgEtaPhi->Fill(deltaPhi, avgEta, weight);
1974 }
1975
1976 // ----- histogram management -----
1977 TH1* 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
1996 TH2* 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
2017 TH3* 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 }