]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskHJetDphi.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
making script executable
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskHJetDphi.cxx
CommitLineData
aac81e6b 1#include <TCanvas.h>
2#include <TChain.h>
3#include <TFormula.h>
4#include <TH1.h>
5#include <TH2.h>
6#include <TH3.h>
7#include <TProfile2D.h>
8#include <THnSparse.h>
9#include <TROOT.h>
10#include <TTree.h>
11#include <TArrayI.h>
12#include <TClonesArray.h>
13#include <TRandom3.h>
14#include <TFile.h>
15#include <TF1.h>
16#include <TLorentzVector.h>
17#include "TKey.h"
18#include "TList.h"
19#include "TSystem.h"
20#include "AliFJWrapper.h"
21#include "AliAODHandler.h"
22#include "AliAODEvent.h"
23#include "AliAODInputHandler.h"
24#include "AliAnalysisManager.h"
25#include "AliAnalysisTask.h"
26#include "AliCentrality.h"
27#include "AliAnalysisTaskHJetDphi.h"
28#include "AliESDEvent.h"
29#include "AliESDInputHandler.h"
30#include "AliESDtrack.h"
31#include "AliESDtrackCuts.h"
32#include "AliVParticle.h"
33#include "AliVTrack.h"
34#include "AliInputEventHandler.h"
35#include "AliMCEvent.h"
36#include "AliStack.h"
37#include "AliGenEventHeader.h"
38#include "AliGenPythiaEventHeader.h"
39#include "AliLog.h"
40#include "AliRhoParameter.h"
41#include "AliAODMCParticle.h"
42#include "AliNamedArrayI.h"
43#include "AliNamedString.h"
44#include "AliPicoTrack.h"
45#include "AliAnalysisTaskFastEmbedding.h"
46#include "AliEmcalJet.h"
47#include "AliAODJet.h"
48#include "AliAODJetEventBackground.h"
abe5a06b 49#include "AliAnalysisHelperJetTasks.h"
aac81e6b 50
51#include <iostream>
52using std::cout;
53using std::cerr;
54using std::endl;
55
56ClassImp(AliAnalysisTaskHJetDphi)
57
58const Double_t pi = TMath::Pi();
aac81e6b 59const Double_t kSector = pi/9;
60
61//________________________________________________________________________
62AliAnalysisTaskHJetDphi::AliAnalysisTaskHJetDphi() :
63 AliAnalysisTaskSE(),
64 fVerbosity(0), fIsEmbedding(kFALSE), fAnaType(0), fPeriod("lhc11h"), fCollisionSystem("PbPb"),
65 fIsMC(kFALSE), fAnalyzeMCTruth(kFALSE), fMC(0),
66 fEvent(0x0), fESD(0x0), fAODIn(0x0), fAODOut(0x0), fAODExtension(0x0),
67 fOfflineTrgMask(AliVEvent::kAny), fTriggerType(-1), fCentrality(-1), fMaxVtxZ(10),
49bf6461 68 fEsdTrkCut(0x0), fEsdHybCut(0x0), fFilterMask(0), fRequireITSRefit(kTRUE), fRequireSharedClsCut(kTRUE),
aac81e6b 69 fIsInit(kFALSE), fNonStdFile(""), fMcParticleArrName(""), fMcParticleArray(0x0), fMcParticlelMap(0x0),
70 fEmbTrkArrName(""), fEmbTrkArray(0x0), fTrackArrName(""), fTrackArray(0x0),
605e0d00 71 fTriggerTrkIndex(-1), fTriggerTrkPt(-1), fSwitchOnAvoidTpcHole(kFALSE), fAvoidTpcHole(0), fCutTPCBoundary(kFALSE), fDistToTPCBoundary(0.),
aac81e6b 72 fMinTrkPt(0.15), fMaxTrkPt(1e4), fMinTrkEta(-0.9), fMaxTrkEta(0.9), fMinTrkPhi(0), fMaxTrkPhi(2*pi),
73 fRadius(0.4), fJetArrName(""), fPLJetArrName(""), fDLJetArrName(""), fJetArray(0x0), fPLJetArray(0x0), fDLJetArray(0x0),
74 fRhoName(""), fRho(0x0), fRhoValue(0), fEvtBkg(0x0), fPtHardBinName(0x0), fPtHardBin(-1),
75 fRunTrkQA(kFALSE), fRunJetQA(kFALSE), fRunSingleInclHJet(kFALSE), fTTtype(0), fTTMinPt(9), fTTMaxPt(10), fJetPtMin(10),
76 fRunPLHJet(kFALSE), fRunDLHJet(kFALSE), fRunLeadTrkQA(kFALSE), fStudyKtEffects(kFALSE), fKtValue(0), fRandom(0),
77 fRunBkgFlow(kFALSE),
78 fOutputList(0x0), fhEventStat(0x0), fhNTrials(0x0), fhPtHardBins(0x0)
79{
80 // Constructor
81
82 // Define input and output slots here
83 // Input slot #0 works with a TChain
84 //DefineInput(0, TChain::Class());
85 //DefineOutput(1, TList::Class());
86 // Output slot #0 id reserved by the base class for AOD
87
88 for(Int_t i=0; i<4; i++)
89 {
90 fhVtxZ[i] = 0x0;
91 fhCentrality[i] = 0x0;
92 fhRhoVsCent[i] = 0x0;
93
94 fhTrkPt[i] = 0x0;
95 fhTrkQA[i] = 0x0;
96 fhTrkPtRes[i] = 0x0;
97 fhTrkPhiRes[i] = 0x0;
98
99 fhNumberOfTT[i] = 0x0;
100 for(Int_t j=0; j<3; j++)
101 {
102 fhJetPt[i][j] = 0x0;
103 fhJetArea[i][j] = 0x0;
104 fhJetQA[i][j] = 0x0;
105
106 fhTTPt[i][j] = 0x0;
107 fHJetPhiCorr[i][j] = 0x0;
108 }
109 fHJetPhiCorrUp[i] = 0x0;
110 fHJetPhiCorrDown[i] = 0x0;
111
112 fhLeadTrkQA[i] = 0x0;
113 fhKtEffects[i] = 0x0;
114 }
115
116 fAODfilterBits[0] = -1;
117 fAODfilterBits[1] = -1;
118}
119
120//________________________________________________________________________
121AliAnalysisTaskHJetDphi::AliAnalysisTaskHJetDphi(const char *name) :
122 AliAnalysisTaskSE(name),
123 fVerbosity(0), fIsEmbedding(kFALSE), fAnaType(0), fPeriod("lhc11h"), fCollisionSystem("PbPb"),
124 fIsMC(kFALSE), fAnalyzeMCTruth(kFALSE), fMC(0),
125 fEvent(0x0), fESD(0x0), fAODIn(0x0), fAODOut(0x0), fAODExtension(0x0),
126 fOfflineTrgMask(AliVEvent::kAny), fTriggerType(-1), fCentrality(-1), fMaxVtxZ(10),
49bf6461 127 fEsdTrkCut(0x0), fEsdHybCut(0x0), fFilterMask(0), fRequireITSRefit(kTRUE), fRequireSharedClsCut(kTRUE),
aac81e6b 128 fIsInit(kFALSE), fNonStdFile(""), fMcParticleArrName(""), fMcParticleArray(0x0), fMcParticlelMap(0x0),
129 fEmbTrkArrName(""), fEmbTrkArray(0x0), fTrackArrName(""), fTrackArray(0x0),
605e0d00 130 fTriggerTrkIndex(-1), fTriggerTrkPt(-1), fSwitchOnAvoidTpcHole(kFALSE), fAvoidTpcHole(0), fCutTPCBoundary(kFALSE), fDistToTPCBoundary(0.),
aac81e6b 131 fMinTrkPt(0.15), fMaxTrkPt(1e4), fMinTrkEta(-0.9), fMaxTrkEta(0.9), fMinTrkPhi(0), fMaxTrkPhi(2*pi),
132 fRadius(0.4), fJetArrName(""), fPLJetArrName(""), fDLJetArrName(""), fJetArray(0x0), fPLJetArray(0x0), fDLJetArray(0x0),
133 fRhoName(""), fRho(0x0), fRhoValue(0), fEvtBkg(0x0), fPtHardBinName(0x0), fPtHardBin(-1),
134 fRunTrkQA(kFALSE), fRunJetQA(kFALSE), fRunSingleInclHJet(kFALSE), fTTtype(0), fTTMinPt(9), fTTMaxPt(10), fJetPtMin(10),
135 fRunPLHJet(kFALSE), fRunDLHJet(kFALSE), fRunLeadTrkQA(kFALSE), fStudyKtEffects(kFALSE), fKtValue(0), fRandom(0),
136 fRunBkgFlow(kFALSE),
137 fOutputList(0x0), fhEventStat(0x0), fhNTrials(0x0), fhPtHardBins(0x0)
138{
139 // Constructor
140
141 // Define input and output slots here
142 // Input slot #0 works with a TChain
143 DefineInput(0, TChain::Class());
144 DefineOutput(1, TList::Class());
145 // Output slot #0 id reserved by the base class for AOD
146
147 for(Int_t i=0; i<4; i++)
148 {
149 fhVtxZ[i] = 0x0;
150 fhCentrality[i] = 0x0;
151 fhRhoVsCent[i] = 0x0;
152
153 fhTrkPt[i] = 0x0;
154 fhTrkQA[i] = 0x0;
155 fhTrkPtRes[i] = 0x0;
156 fhTrkPhiRes[i] = 0x0;
157
158 fhNumberOfTT[i] = 0x0;
159 for(Int_t j=0; j<3; j++)
160 {
161 fhJetPt[i][j] = 0x0;
162 fhJetArea[i][j] = 0x0;
163 fhJetQA[i][j] = 0x0;
164
165 fhTTPt[i][j] = 0x0;
166 fHJetPhiCorr[i][j] = 0x0;
167 }
168
169 fHJetPhiCorrUp[i] = 0x0;
170 fHJetPhiCorrDown[i] = 0x0;
171
172 fhLeadTrkQA[i] = 0x0;
173 fhKtEffects[i] = 0x0;
174 }
175
176 fAODfilterBits[0] = -1;
177 fAODfilterBits[1] = -1;
178}
179
180
181//________________________________________________________________________
182AliAnalysisTaskHJetDphi::~AliAnalysisTaskHJetDphi()
183{
184 //Destructor
185 if(fRandom) delete fRandom;
186 if(fEsdTrkCut) delete fEsdTrkCut;
187 if(fEsdHybCut) delete fEsdHybCut;
188 if(fOutputList) { fOutputList->Delete(); delete fOutputList;}
189}
190
191//________________________________________________________________________
192void AliAnalysisTaskHJetDphi::UserCreateOutputObjects()
193{
194 // Create histograms
195
196 const Int_t nTrkPtBins = 100;
197 const Float_t lowTrkPtBin=0, upTrkPtBin=100;
198 const Int_t nJetPtBins = 300;
199 const Float_t lowJetPtBin=-100, upJetPtBin=200;
200
201 // track QA
202 const Int_t dimTrkqa = 4;
203 const Int_t nBinsTrkqa[dimTrkqa] = {nTrkPtBins/5, 36, 40, 10};
204 const Double_t lowBinTrkqa[dimTrkqa] = {lowTrkPtBin, 0, -1, 0};
205 const Double_t hiBinTrkqa[dimTrkqa] = {upTrkPtBin, 360, 1, 10};
206
207 const Int_t dimTrkRes = 5;
208 const Int_t nBinsTrkRes[dimTrkRes] = {nTrkPtBins, 50, 50, 3, 10};
209 const Double_t lowBinTrkRes[dimTrkRes] = {lowTrkPtBin, 0, 0, 0, 0};
210 const Double_t hiBinTrkRes[dimTrkRes] = {upTrkPtBin, 0.5, 0.5, 3, 10};
211
212 const Int_t dimPhiRes = 4;
213 const Int_t nBinsPhiRes[dimPhiRes] = {nTrkPtBins, 200, 3, 10};
214 const Double_t lowBinPhiRes[dimPhiRes] = {lowTrkPtBin, -0.00995, 0, 0};
215 const Double_t hiBinPhiRes[dimPhiRes] = {upTrkPtBin, 0.01005, 3, 10};
216
217 // jet QA
218 const Int_t dimJetpt = 4;
219 const Int_t nBinsJetpt[dimJetpt] = {nJetPtBins, 300, 10, 10};
220 const Double_t lowBinJetpt[dimJetpt] = {lowJetPtBin, 0, 0, 0};
221 const Double_t hiBinJetpt[dimJetpt] = {upJetPtBin, 300, 10, 10};
222
223 const Int_t dimJetA = 4;
224 const Int_t nBinsJetA[dimJetA] = {nJetPtBins, 100, 10, 10};
225 const Double_t lowBinJetA[dimJetA] = {lowJetPtBin, 0, 0, 0};
226 const Double_t hiBinJetA[dimJetA] = {upJetPtBin, 1, 10, 10};
227
228 const Int_t dimJetqa = 7;
605e0d00 229 const Int_t nBinsJetqa[dimJetqa] = {nJetPtBins/5, 36, 24, 6, 100, 10, 11};
aac81e6b 230 const Double_t lowBinJetqa[dimJetqa] = {lowJetPtBin, 0, -0.6, 0, 0, 0, 0};
231 const Double_t hiBinJetqa[dimJetqa] = {upJetPtBin, 360, 0.6, 1.2, 500, 10, 11};
232
233 // h-jet analysis
234 const Int_t dimTT = 3;
235 const Int_t nBinsTT[dimTT] = {nTrkPtBins, 10, 11};
236 const Double_t lowBinTT[dimTT] = {lowTrkPtBin, 0, 0};
237 const Double_t hiBinTT[dimTT] = {upTrkPtBin, 100, 11};
238
605e0d00 239 const Int_t dimCor = 8;
240 const Int_t nBinsCor[dimCor] = {nTrkPtBins, nJetPtBins, 140, 6, 10, 40, 11, 10};
241 const Double_t lowBinCor[dimCor] = {lowTrkPtBin,lowJetPtBin, pi-4.95, 0, 0, -1.95, 0, 0};
242 const Double_t hiBinCor[dimCor] = {upTrkPtBin, upJetPtBin, pi+2.05, 1.2, 100, 2.05, 11, 10};
aac81e6b 243
244 // Leading track QA
245 const Int_t dimLeadTrkqa = 5;
246 const Int_t nBinsLeadTrkqa[dimLeadTrkqa] = {nTrkPtBins, 200, 80, 55, 10};
247 const Double_t lowBinLeadTrkqa[dimLeadTrkqa] = {lowTrkPtBin, 0, -0.4, 0, 0};
248 const Double_t hiBinLeadTrkqa[dimLeadTrkqa] = {upTrkPtBin, 200, 0.4, 1.1, 100};
249
250 // kt effects
251 const Int_t dimKt = 5;
252 const Int_t nBinsKt[dimKt] = {nTrkPtBins, 20, 81, 10, 11};
253 const Double_t lowBinKt[dimKt] = {lowTrkPtBin, 0, -40.5, 0, 0};
254 const Double_t hiBinKt[dimKt] = {upTrkPtBin, 200, 40.5, 10, 11};
255
aac81e6b 256 // Called once
257 OpenFile(1);
258 if(!fOutputList) fOutputList = new TList;
259 fOutputList->SetOwner(kTRUE);
260
261 Bool_t oldStatus = TH1::AddDirectoryStatus();
262 if(fAnaType==1) TH1::AddDirectory(kFALSE);
263
abe5a06b 264 fhEventStat = new TH1F("fhEventStat","Event statistics for jet analysis",8,0,8);
aac81e6b 265 fhEventStat->GetXaxis()->SetBinLabel(1,"All");
266 fhEventStat->GetXaxis()->SetBinLabel(2,"PS");
267 fhEventStat->GetXaxis()->SetBinLabel(3,"Vtx");
268 fhEventStat->GetXaxis()->SetBinLabel(4,"Vtx+10cm");
269 fhEventStat->GetXaxis()->SetBinLabel(5,"kMB");
270 fhEventStat->GetXaxis()->SetBinLabel(6,"kCentral");
271 fhEventStat->GetXaxis()->SetBinLabel(7,"kSemiCentral");
abe5a06b 272 fhEventStat->GetXaxis()->SetBinLabel(8,"kJetService");
aac81e6b 273 fOutputList->Add(fhEventStat);
274
275 fhNTrials = new TH1F("fhNTrials","Number of trials",1,0,1);
276 fOutputList->Add(fhNTrials);
277
278 fhPtHardBins = new TH1F("fhPtHardBins","Number of events in each pT hard bin",11,0,11);
279 fOutputList->Add(fhPtHardBins);
280
281 const char *triggerName[4] = {"kMB","kEGA","kEJE","MC"};
282 const char *dataType[3] = {"", "_DL","_PL"};
283 const char *dataName[3] = {"Data","DL","PL"};
284
285 Double_t newbins[7] = {0,0.07,0.2,0.4,0.6,0.8,1};
286
287 for(Int_t i=0; i<4; i++)
288 {
289 if( fAnalyzeMCTruth )
290 {
291 if(i!=3) continue;
292 }
293 else
294 {
295 if( fPeriod=="lhc11a" && i>1 ) continue;
296 if( fPeriod=="lhc10h" && i!=0 ) continue;
297 if( fPeriod=="lhc11h" && i!=0 ) continue;
298 if( fPeriod.Contains("lhc12a15a") && i!=0 ) continue;
299 if( fPeriod.Contains("lhc12a15e") && i!=0 ) continue;
300 }
301
302 fhVtxZ[i] = new TH1F(Form("%s_fhVtxZ",triggerName[i]),Form("%s: z distribution of event vertexz;z(cm)",triggerName[i]),400,-20,20);
303 fOutputList->Add(fhVtxZ[i]);
304
305 fhCentrality[i] = new TH1F(Form("%s_fhCentrality",triggerName[i]),Form("%s: Event centrality;centrality",triggerName[i]),100,0,100);
306 fOutputList->Add(fhCentrality[i]);
307
308 fhRhoVsCent[i] = new TH2F(Form("%s_fhRhoVsCent",triggerName[i]),Form("%s: Rho vs centrality (R=%1.1f);centrality;Rho",triggerName[i],fRadius),100,0,100,300,0,300);
309 fOutputList->Add(fhRhoVsCent[i]);
310
311 if(fRunTrkQA)
312 {
605e0d00 313 fhTrkPt[i] = new TH2F(Form("%s_fhTrkPt",triggerName[i]),Form("%s: Track p_{T} vs centrality;p_{T}^{track} (GeV/c);Centrality",triggerName[i]),nTrkPtBins,lowTrkPtBin,upTrkPtBin,10,0,100);
aac81e6b 314 fOutputList->Add(fhTrkPt[i]);
315
316 fhTrkQA[i] = new THnSparseF(Form("%s_fhTrkQA",triggerName[i]),Form("%s: track p_{T} vs #phi vs #eta vs centrality;p_{T,track} (GeV/c);#phi;#eta;centrality",triggerName[i]),dimTrkqa,nBinsTrkqa,lowBinTrkqa,hiBinTrkqa);
317 fOutputList->Add(fhTrkQA[i]);
318
319 fhTrkPtRes[i] = new THnSparseF(Form("%s_fhTrkPtRes",triggerName[i]),Form("%s: track p_{T} vs resolution vs (p_{T}^{gen}-p_{T}^{rec})/p_{T}^{gen} vs type vs centrality;p_{T,track} (GeV/c);#sigma(p_{T})/p_{T};type;centrality",triggerName[i]),dimTrkRes,nBinsTrkRes,lowBinTrkRes,hiBinTrkRes);
320 fOutputList->Add(fhTrkPtRes[i]);
321
322 fhTrkPhiRes[i] = new THnSparseF(Form("%s_fhTrkPhiRes",triggerName[i]),Form("%s: track p_{T} vs #varphi^{gen}-#varphi^{rec} vs type vs centrality;p_{T,track} (GeV/c);#Delta#varphi;type;centrality",triggerName[i]),dimPhiRes,nBinsPhiRes,lowBinPhiRes,hiBinPhiRes);
323 fOutputList->Add(fhTrkPhiRes[i]);
324 }
325
326 for(Int_t j=0; j<3; j++)
327 {
328 if(!fRunDLHJet && j==1) continue;
329 if(!fRunPLHJet && j==2) continue;
330 if(fRunJetQA)
331 {
332 fhJetPt[i][j] = new THnSparseF(Form("%s_fhJetPt%s",triggerName[i],dataType[j]),Form("%s-%s: jet p_{T} vs raw jet p_{T} vs centrality vs pt hard bin (R=%1.1f);p_{T,jet}^{ch} (GeV/c);p_{T,jet}^{raw} (GeV/c);centrality;ptHardBin",dataName[j],triggerName[i],fRadius),dimJetpt,nBinsJetpt,lowBinJetpt,hiBinJetpt);
333 fOutputList->Add(fhJetPt[i][j]);
334
335 fhJetArea[i][j] = new THnSparseF(Form("%s_fhJetArea%s",triggerName[i],dataType[j]),Form("%s-%s: jet p_{T} vs area vs centrality vs pt hard bin (R=%1.1f);p_{T,jet}^{ch} (GeV/c);area;centrality;ptHardBin",dataName[j],triggerName[i],fRadius),dimJetA,nBinsJetA,lowBinJetA,hiBinJetA);
336 fOutputList->Add(fhJetArea[i][j]);
337
338 fhJetQA[i][j] = new THnSparseF(Form("%s_fhJetQA%s",triggerName[i],dataType[j]),Form("%s-%s: jet p_{T} vs #phi vs #eta vs area vs # of constituents vs centrality vs pt hard bin (R=%1.1f);p_{T,jet}^{ch} (GeV/c);#phi;#eta;area;# of constituents;centrality;ptHardBin",dataName[j],triggerName[i],fRadius),dimJetqa,nBinsJetqa,lowBinJetqa,hiBinJetqa);
339 fhJetQA[i][j]->SetBinEdges(3,newbins);
340 fOutputList->Add(fhJetQA[i][j]);
341 }
342
343 if(fRunSingleInclHJet)
344 {
345 fhTTPt[i][j] = new THnSparseF(Form("%s_fhTTPt%s",triggerName[i],dataType[j]),Form("%s-%s: TT p_{T} vs centrality vs pT hard bin;p_{T,TT}^{ch} (GeV/c);centrality;ptHardBin",dataName[j],triggerName[i]),dimTT,nBinsTT,lowBinTT,hiBinTT);
346 fOutputList->Add(fhTTPt[i][j]);
347
605e0d00 348 fHJetPhiCorr[i][j] = new THnSparseF(Form("%s_fHJetPhiCorr%s",triggerName[i],dataType[j]),Form("%s-%s: single inclusive TT p_{T} vs recoil jet p_{T} vs #Delta#varphi vs area vs centrality vs #Delta#eta vs pt hard bin vs event bin (R=%1.1f);TT p_{T} (GeV/c);p_{T,jet}^{ch} (GeV/c);#Delta#varphi;area;centrality;#Delta#eta;ptHardBin;EventBin",dataName[j],triggerName[i],fRadius),dimCor,nBinsCor,lowBinCor,hiBinCor);
aac81e6b 349 fHJetPhiCorr[i][j]->SetBinEdges(3,newbins);
350 fOutputList->Add(fHJetPhiCorr[i][j]);
351 }
352 }
353
354
355 if(fRunSingleInclHJet)
356 {
357 fhNumberOfTT[i] = new TH1F(Form("%s_fhNumberOfTT",triggerName[i]), Form("%s: number of TT",triggerName[i]),6,0,6);
358 fOutputList->Add(fhNumberOfTT[i]);
359
aac81e6b 360
0330533a 361 if(fRunBkgFlow)
362 {
605e0d00 363 fHJetPhiCorrUp[i] = new THnSparseF(Form("%s_fHJetPhiCorrUp",triggerName[i]),Form("%s: single inclusive TT p_{T} vs recoil jet p_{T} vs #Delta#varphi vs area vs centrality vs #Delta#eta vs pt hard bin vs event bin (R=%1.1f);TT p_{T} (GeV/c);p_{T,jet}^{ch} (GeV/c);#Delta#varphi;area;centrality;#Delta#eta;ptHardBin;EventBin",triggerName[i],fRadius),dimCor,nBinsCor,lowBinCor,hiBinCor);
0330533a 364 fOutputList->Add(fHJetPhiCorrUp[i]);
365
605e0d00 366 fHJetPhiCorrDown[i] = new THnSparseF(Form("%s_fHJetPhiCorrDown",triggerName[i]),Form("%s: single inclusive TT p_{T} vs recoil jet p_{T} vs #Delta#varphi vs area vs centrality vs #Delta#eta vs pt hard bin vs event bin (R=%1.1f);TT p_{T} (GeV/c);p_{T,jet}^{ch} (GeV/c);#Delta#varphi;area;centrality;#Delta#eta;ptHardBin;EventBin",triggerName[i],fRadius),dimCor,nBinsCor,lowBinCor,hiBinCor);
0330533a 367 fOutputList->Add(fHJetPhiCorrDown[i]);
368 }
aac81e6b 369 }
370
371 if(fRunLeadTrkQA)
372 {
373 fhLeadTrkQA[i] = new THnSparseF(Form("%s_fhLeadTrkQA",triggerName[i]),Form("%s: p_{T,trk}^{leading} vs p_{T,jet}^{full} vs #Delta#varphi vs z vs centrality;p_{T,trk}^{leading} (GeV/c); p_{T,jet}^{full} (GeV/c);#Delta#varphi;z;centrality",triggerName[i]),dimLeadTrkqa,nBinsLeadTrkqa,lowBinLeadTrkqa,hiBinLeadTrkqa);
374 fOutputList->Add(fhLeadTrkQA[i]);
375 }
376
377 if(fStudyKtEffects)
378 {
379 fhKtEffects[i] = new THnSparseF(Form("%s_fhKtEffects",triggerName[i]),Form("%s: TT p_{T} vs recoil jet p_{T} vs k_{t} vs centrality vs pt hard bin (R=%1.1f);TT p_{T} (GeV/c);p_{T,jet}^{ch} (GeV/c);k_{t} (GeV/c);centrality;ptHardBin",triggerName[i],fRadius),dimKt,nBinsKt,lowBinKt,hiBinKt);
380 fOutputList->Add(fhKtEffects[i]);
381 }
382 }
383
384 //error calculation in THnSparse
385 Int_t nObj = fOutputList->GetEntries();
386 for(Int_t i=0; i<nObj; i++)
387 {
388 TObject *obj = (TObject*) fOutputList->At(i);
389 if (obj->IsA()->InheritsFrom( "THnSparse" ))
390 {
391 THnSparseF *hn = (THnSparseF*)obj;
392 hn->Sumw2();
393 }
394 }
395
396 if(fRunTrkQA)
397 {
398 fEsdTrkCut = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE,1);
399 fEsdTrkCut->SetMaxDCAToVertexXY(2.4);
400 fEsdTrkCut->SetMaxDCAToVertexZ(3.2);
401 fEsdTrkCut->SetDCAToVertex2D(kTRUE);
402 fEsdTrkCut->SetMaxChi2TPCConstrainedGlobal(36);
403 fEsdTrkCut->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
404
405 fEsdHybCut = new AliESDtrackCuts(*fEsdTrkCut);
406 fEsdHybCut->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kNone);
407 }
408
409 fRandom = new TRandom3(0);
410
411 PrintConfig();
412
413 if(fAnaType==1) TH1::AddDirectory(oldStatus);
414 PostData(1, fOutputList);
415}
416
417//________________________________________________________________________
418Bool_t AliAnalysisTaskHJetDphi::UserNotify()
419{
420 AliInfo("User Nofity");
421
422 Int_t runNumber = InputEvent()->GetRunNumber();
605e0d00 423
424 fAvoidTpcHole = 0;
aac81e6b 425 if(fSwitchOnAvoidTpcHole)
426 {
605e0d00 427 Int_t runs_iroc[28] = {169975, 169981, 170038, 170040, 170083, 170084, 170085, 170088, 170089, 170091, 170152, 170155, 170159, 170163, 170193, 170195, 170203, 170204, 170205, 170228, 170230, 170264, 170268, 170269, 170270, 170306, 170308, 170309};
428 Int_t runs_oroc[23] = {169591, 169590, 169588, 169587, 169586, 169584, 169557, 169555, 169554, 169553, 169550, 169515, 169512, 169506, 169504, 169498, 169475, 169420, 169418, 169099, 169040, 169045, 169044};
429
aac81e6b 430 for(Int_t i=0; i<28; i++)
431 {
605e0d00 432 if(runNumber==runs_iroc[i])
433 {
434 fAvoidTpcHole = 1;
435 break;
436 }
437 }
438 for(Int_t i=0; i<23; i++)
439 {
440 if(runNumber==runs_oroc[i])
aac81e6b 441 {
605e0d00 442 fAvoidTpcHole = 2;
aac81e6b 443 break;
444 }
445 }
aac81e6b 446 }
447
448 if(fIsMC)
449 {
450 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
451 TFile *currFile = tree->GetCurrentFile();
452 TString fileName(currFile->GetName());
453 if(fileName.Contains("root_archive.zip#"))
454 {
455 Ssiz_t pos = fileName.Index("#",0,TString::kExact);
456 fileName.Replace(pos+1,20,"");
457 }
458 else
459 {
460 fileName.ReplaceAll(gSystem->BaseName(fileName.Data()),"");
461 }
462
463 TFile *fxsec = TFile::Open(Form("%s%s",fileName.Data(),"pyxsec_hists.root"),"read");
464 if(fxsec)
465 {
466 TKey *key = (TKey*)fxsec->GetListOfKeys()->At(0);
467 if(key)
468 {
469 TList *list = dynamic_cast<TList*>(key->ReadObj());
470 if(list)
471 {
472 fhNTrials->Fill(0.5, ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1));
473 }
474 else
475 return kFALSE;
476 }
477 else
478 return kFALSE;
479 }
480 else
481 {
482 fxsec = TFile::Open(Form("%s%s",fileName.Data(),"pyxsec.root"),"read");
483 TTree *xtree = (TTree*)fxsec->Get("Xsection");
484 if(xtree)
485 {
486 UInt_t ntrials = 0;
487 xtree->SetBranchAddress("ntrials",&ntrials);
488 xtree->GetEntry(0);
489 fhNTrials->Fill(0.5, ntrials);
490 }
491 else
492 return kFALSE;
493 }
494 fxsec->Close();
495 }
496 return kTRUE;
497}
498
499//________________________________________________________________________
500void AliAnalysisTaskHJetDphi::UserExec(Option_t *)
501{
502 // Main loop, called for each event.
503
504 fTriggerType = -1;
505
506 // Get pointers to input events
507 fEvent = InputEvent();
508 if (!fEvent)
509 {
510 AliError("Input event not available");
511 return;
512 }
513
514 if(fIsMC)
515 {
516 fMC = MCEvent();
517 if (!fMC)
518 {
519 AliError("MC event available");
520 return;
521 }
522 }
523
524 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
525 if (!fESD)
526 {
aac81e6b 527 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
528 }
0330533a 529 else
530 {
531 if(fAnaType==1) fEvent = AODEvent();
532 }
aac81e6b 533 if(fAnaType==1)
534 {
535 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
536 if(fNonStdFile.Length()!=0)
537 {
538 // case that we have an AOD extension we need can fetch the jets from the extended output
539 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
540 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
541 if(!fAODExtension)
542 {
543 if(fVerbosity>1) printf("AODExtension not found for %s",fNonStdFile.Data());
544 }
545 }
546 }
547
548 fhEventStat->Fill(0.5);
549
0330533a 550 if(fVerbosity>1)
551 {
49bf6461 552 TList *list = 0x0;
553 if(fAnaType==0) list = fEvent->GetList();
554 else list = fAODOut->GetList();
0330533a 555 for(Int_t i=0; i<list->GetEntries(); i++)
556 {
557 TObject *obj = (TObject*)list->At(i);
558 cout<<i<<": "<<obj->GetName()<<" : "<<obj->ClassName()<<endl;
559 }
560 }
561 // Retrieve arraies from memory
562 if(!fIsInit) fIsInit = RetrieveArraies();
563 if(!fIsInit) return;
564
aac81e6b 565 if(fIsEmbedding)
566 {
567 if(fAnaType==0)
568 {
569 TString fileName = fPtHardBinName->GetString();
570 fileName.Remove(0,51);
571 fileName.Remove(fileName.Index("/"));
572 fPtHardBin = fileName.Atoi();
573 }
574 if(fAnaType==1)
575 {
576 Double_t pthard = AliAnalysisTaskFastEmbedding::GetPtHard();
577 fPtHardBin = GetPtHardBin(pthard);
578 }
579 }
580
581 // physics selection
582 UInt_t trigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
583 if(fOfflineTrgMask==AliVEvent::kAny) fTriggerType = 0;
584 else
585 {
586 if(trigger & fOfflineTrgMask) fTriggerType=0;
587
588 if(fPeriod.Contains("lhc11a",TString::kIgnoreCase))
589 {
590 if (trigger & AliVEvent::kMB) { fTriggerType=0; }
591 if (trigger & AliVEvent::kEMC1) { fTriggerType=1; }
592 }
593 }
594 if(fTriggerType==-1) return;
595 if(fAnalyzeMCTruth) fTriggerType = 3;
596 fhEventStat->Fill(1.5);
597
598 // Vertex cut
599 const AliVVertex* vtx = fEvent->GetPrimaryVertex();
600 if (!vtx || vtx->GetNContributors()<1) return;
601 fhEventStat->Fill(2.5);
602 fhVtxZ[fTriggerType]->Fill(vtx->GetZ());
603 if (TMath::Abs(vtx->GetZ())>fMaxVtxZ) return;
604 fhEventStat->Fill(3.5);
605
606 if(fTriggerType==0)
607 {
608 if(trigger & AliVEvent::kCentral) fhEventStat->Fill(5.5);
609 else if (trigger & AliVEvent::kCentral) fhEventStat->Fill(6.5);
610 else fhEventStat->Fill(4.5);
611 }
abe5a06b
CL		 612
613 if(!AliAnalysisHelperJetTasks::Selected()) return;
614 fhEventStat->Fill(7.5);
aac81e6b 615
616 // GetCentrality
617 if(fCollisionSystem=="PbPb")
618 {
619 AliCentrality *centrality = fEvent->GetCentrality();
620 if (centrality) fCentrality = centrality->GetCentralityPercentile("V0M");
621 else fCentrality = 99;
622 }
623 else if(fCollisionSystem=="pp")
624 fCentrality = 0;
625 fhCentrality[fTriggerType]->Fill(fCentrality);
626
aac81e6b 627 // Get Rho value
628 if(fCollisionSystem=="PbPb")
629 {
630 if(fAnaType==0) fRhoValue = fRho->GetVal();
631 if(fAnaType==1) fRhoValue = fEvtBkg->GetBackground(0);
632 }
633 else if(fCollisionSystem=="pp")
634 {
635 fRhoValue = 0;
636 }
637 fhRhoVsCent[fTriggerType]->Fill(fCentrality,fRhoValue);
638 fhPtHardBins->Fill(fPtHardBin);
639
640 if(fRunSingleInclHJet)
641 {
642 Double_t trigPt = -1, trigPhi = -999, trigEta = -999;
643 Int_t nTrig = FindSingleIncTrigger(fTrackArray, trigPt, trigPhi, trigEta, fIsEmbedding);
644 if(nTrig>0) fhNumberOfTT[fTriggerType]->Fill(nTrig);
645 RunSingleInclHJetCorr(trigPt, trigPhi, trigEta, fJetArray, fRhoValue, fhTTPt[fTriggerType][0], fHJetPhiCorr[fTriggerType][0]);
646 if(fRunBkgFlow)
647 {
648 RunSingleInclHJetCorr(trigPt, trigPhi, trigEta, fJetArray, fRhoValue+1, 0x0, fHJetPhiCorrUp[fTriggerType]);
649 RunSingleInclHJetCorr(trigPt, trigPhi, trigEta, fJetArray, fRhoValue-1, 0x0, fHJetPhiCorrDown[fTriggerType]);
650 }
651
652 if(fIsEmbedding)
653 {
654 if(fRunDLHJet)
655 {
656 FindSingleIncTrigger(fTrackArray, trigPt, trigPhi, trigEta, 1);
657 RunSingleInclHJetCorr(trigPt, trigPhi, trigEta, fDLJetArray, 0, fhTTPt[fTriggerType][1], fHJetPhiCorr[fTriggerType][1]);
658 }
659
660 if(fRunPLHJet)
661 {
662 FindSingleIncTrigger(fMcParticleArray, trigPt, trigPhi, trigEta, 2);
663 RunSingleInclHJetCorr(trigPt, trigPhi, trigEta, fPLJetArray, 0, fhTTPt[fTriggerType][2], fHJetPhiCorr[fTriggerType][2]);
664 }
665 }
666 }
667
668 if(fRunJetQA)
669 {
670 RunJetQA(fJetArray, fRhoValue, fhJetPt[fTriggerType][0], fhJetArea[fTriggerType][0], fhJetQA[fTriggerType][0]);
671 if(fIsEmbedding)
672 {
673 if(fRunDLHJet) RunJetQA(fDLJetArray, 0, fhJetPt[fTriggerType][1], fhJetArea[fTriggerType][1], fhJetQA[fTriggerType][1]);
674 if(fRunPLHJet) RunJetQA(fPLJetArray, 0, fhJetPt[fTriggerType][2], fhJetArea[fTriggerType][2], fhJetQA[fTriggerType][2]);
675 }
676 }
677
678 if(!fIsEmbedding)
679 {
680 if(fRunTrkQA) RunTrackQA();
681 if(fRunLeadTrkQA) RunLeadTrkQA();
682 }
683
684 if(fStudyKtEffects) StudyKtEffects();
685
686 PostData(1, fOutputList);
687 return;
688}
689
690//________________________________________________________________________
691Int_t AliAnalysisTaskHJetDphi::FindSingleIncTrigger(const TClonesArray *trackArray, Double_t &trigPt, Double_t &trigPhi, Double_t &trigEta, const Int_t arrayType)
692{
693 Int_t trigIndex = -1;
694 trigPt = -1;
695 trigPhi = -999;
696 trigEta = -999;
697
698 // arrayType: 0 -> data, 1 -> embedding, 2 -> MC
699 if(!trackArray) return 0;
700
701 Int_t nTT = 0, counter = 0;
702 Int_t ntracks = trackArray->GetEntries();
703 TArrayI arr;
704 arr.Set(ntracks);
705 for(Int_t it=0; it<ntracks; it++)
706 {
707 AliVParticle *t = dynamic_cast<AliVParticle*>(trackArray->At(it));
708 if(!t || t->Charge()==0 || !AcceptTrack(t)) continue;
709
710 if(fAnaType==0 && arrayType==1 && t->GetLabel()==0) continue;
711
712 if(fAnaType==1 && arrayType<2 && !IsGoodAODtrack(t) ) continue;
713
714 Double_t pt = t->Pt();
715 if(fTTtype==0)
716 {
717 if (pt<fTTMaxPt && pt>=fTTMinPt)
718 {
719 nTT++;
720 arr.AddAt(it,counter);
721 counter++;
722 }
723 }
724 }
725 arr.Set(counter);
726
727 if(fTTtype==0)
728 {
729 if(counter==0) trigIndex = -1;
730 else if(counter==1) trigIndex = arr.At(0);
731 else
732 {
733 fRandom->SetSeed(arr.At(0)); //make this random selection reproducible
734 Double_t pro = fRandom->Uniform() * counter;
735 trigIndex = arr.At(TMath::FloorNint(pro));
736 }
737 arr.Reset();
738 }
739
0330533a 740 if(trigIndex>-1)
aac81e6b 741 {
742 AliVParticle *tt = (AliVParticle*) fTrackArray->At(trigIndex);
743 trigPt = tt->Pt();
744 trigPhi = tt->Phi();
745 trigEta = tt->Eta();
605e0d00 746
747 if(fSwitchOnAvoidTpcHole)
748 {
749 if(fAvoidTpcHole==1 && !(trigPhi>3.89 && trigPhi<5.53)) trigIndex = -1;
750 if(fAvoidTpcHole==2 && !(trigPhi>2.45 && trigPhi<3.44)) trigIndex = -1;
751 }
aac81e6b 752
753 if(fCutTPCBoundary)
754 {
755 Double_t phiDist = trigPhi - TMath::FloorNint(trigPhi/kSector)*kSector;
756 if(phiDist<fDistToTPCBoundary || phiDist>kSector-fDistToTPCBoundary)
757 {
758 trigIndex = -1;
759 }
760 }
761 }
0330533a 762
aac81e6b 763 if(trigIndex==-1) { trigPt = -1; trigPhi = -999; trigEta = -999;}
0330533a 764 if(arrayType<2) fTriggerTrkIndex = trigIndex;
aac81e6b 765 return nTT;
766}
767
768//________________________________________________________________________
769void AliAnalysisTaskHJetDphi::RunSingleInclHJetCorr(Double_t trigPt, Double_t trigPhi, Double_t trigEta, const TClonesArray *jetArray, Double_t rho, THnSparse *hTT, THnSparse *hn)
770{
771 if(trigPt<0 || !fJetArray) return;
772
773 if(hTT)
774 {
0330533a 775 Double_t fillTT[] = {trigPt, fCentrality, (Double_t)fPtHardBin};
aac81e6b 776 hTT->Fill(fillTT);
777 }
778
779 Int_t nJets = jetArray->GetEntries();
0330533a 780 Double_t jetPt = 0, jetEta = 0, jetPhi = 0, jetArea = 0;
aac81e6b 781 for(Int_t ij=0; ij<nJets; ij++)
782 {
783 if(fAnaType==0)
784 {
785 AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(jetArray->At(ij));
786 jetPt = jet->Pt();
787 jetEta = jet->Eta();
788 jetPhi = jet->Phi();
789 jetArea = jet->Area();
790 }
0330533a 791 else if(fAnaType==1)
aac81e6b 792 {
793 AliAODJet* jet = dynamic_cast<AliAODJet*>(jetArray->At(ij));
794 jetPt = jet->Pt();
795 jetEta = jet->Eta();
796 jetPhi = jet->Phi();
797 jetArea = jet->EffectiveAreaCharged();
798 }
0330533a 799 else
800 return;
aac81e6b 801 if(!IsGoodJet(jetEta)) continue; // eta cut
802 Double_t dPhi = CalculateDPhi(trigPhi,jetPhi);
803
804 Double_t jetPtCorr = jetPt-rho*jetArea;
805 if(jetPtCorr>fJetPtMin)
806 {
c7ac2d0e 807 Double_t fill[] = {trigPt,jetPtCorr,dPhi,jetArea,fCentrality,trigEta-jetEta, (Double_t)fPtHardBin,static_cast<Double_t>(Entry()%10)};
aac81e6b 808 hn->Fill(fill);
809 }
810 }
811}
812
813
814//________________________________________________________________________
815void AliAnalysisTaskHJetDphi::RunTrackQA()
816{
817 if(fIsEmbedding) return;
818 if(!fTrackArray) return;
819
820 Int_t ntracks = fTrackArray->GetEntries();
821 for(Int_t it=0; it<ntracks; it++)
822 {
823 AliVParticle *t = dynamic_cast<AliVParticle*>(fTrackArray->At(it));
824 if(!t || t->Charge()==0 || !AcceptTrack(t)) continue;
825 if(fAnaType==1 && !IsGoodAODtrack(t)) continue;
826 fhTrkPt[fTriggerType]->Fill(t->Pt(),fCentrality);
827 Double_t phi = t->Phi();
828 if(phi<0) phi += 2*pi;
829 Double_t fill[] = {t->Pt(), phi*TMath::RadToDeg(), t->Eta(), fCentrality};
830 fhTrkQA[fTriggerType]->Fill(fill);
831 }
832
833 if(fESD)
834 {
835 Int_t ntrack = fESD->GetNumberOfTracks();
836 for(Int_t itr=0; itr<ntrack; itr++)
837 {
838 AliESDtrack *esdtrack = fESD->GetTrack(itr);
839 if(!esdtrack || TMath::Abs(esdtrack->Eta())>0.9 )continue;
840 Int_t type = -1;
841 if(fEsdTrkCut->AcceptTrack(esdtrack))
842 type = 0;
843 else if(fEsdHybCut->AcceptTrack(esdtrack))
844 type = 1;
845 else
846 continue;
847
848 Double_t resolution = -1;
849 Int_t label = esdtrack->GetLabel();
850 if(label>0 && fMC)
851 {
852 AliMCParticle *part = (AliMCParticle*)fMC->GetTrack(label);
853 TParticle *tpart = (TParticle*)part->Particle();
854 if(TMath::Abs(tpart->GetPdgCode())==211)
855 {
0330533a 856 Double_t fillPhiRes[] = {esdtrack->Pt(),part->Phi()-esdtrack->Phi(),(Double_t)type,fCentrality};
aac81e6b 857 fhTrkPhiRes[fTriggerType]->Fill(fillPhiRes);
858 }
859 resolution = (part->Pt()-esdtrack->Pt())/part->Pt();
860 }
0330533a 861 Double_t fillPtRes[] = {esdtrack->Pt(),esdtrack->Pt()*TMath::Sqrt(esdtrack->GetSigma1Pt2()),resolution,(Double_t)type,fCentrality};
aac81e6b 862 fhTrkPtRes[fTriggerType]->Fill(fillPtRes);
863 }
864 }
865 else if(fAODIn)
866 {
867 ntracks = fAODIn->GetNumberOfTracks();
868 Int_t type = -1;
869 for(Int_t itrack=0; itrack<ntracks; itrack++)
870 {
871 AliAODTrack *aodtrack = (AliAODTrack*)fAODIn->GetTrack(itrack);
872 if(!aodtrack || !AcceptTrack(dynamic_cast<AliVParticle*>(aodtrack)) ) continue;
873 if (fAODfilterBits[0] < 0)
874 {
875 if (aodtrack->IsHybridGlobalConstrainedGlobal())
876 type = 3;
877 else
878 continue;
879 }
880 else
881 {
882 if (aodtrack->TestFilterBit(fAODfilterBits[0])) type = 0;
883 else if (aodtrack->TestFilterBit(fAODfilterBits[1]) && (aodtrack->GetStatus()&AliESDtrack::kITSrefit)!=0) type = 1;
884 else continue;
885 }
886 Double_t sigma1Pt2 = GetAODTrackPtRes(aodtrack);
887 Double_t resolution = -1;
888 Int_t label = aodtrack->GetLabel();
889 if(label>0 && fMC)
890 {
891 AliMCParticle *part = (AliMCParticle*)fMC->GetTrack(label);
892 resolution = (part->Pt()-aodtrack->Pt())/part->Pt();
0330533a 893 Double_t fillPhiRes[] = {aodtrack->Pt(),part->Phi()-aodtrack->Phi(),(Double_t)type,fCentrality};
aac81e6b 894 fhTrkPhiRes[fTriggerType]->Fill(fillPhiRes);
895 }
896 if(sigma1Pt2>0)
897 {
0330533a 898 Double_t fillPtRes[5] = {aodtrack->Pt(),aodtrack->Pt()*TMath::Sqrt(sigma1Pt2),resolution,(Double_t)type,fCentrality};
aac81e6b 899 fhTrkPtRes[fTriggerType]->Fill(fillPtRes);
900 }
901 }
902 }
903}
904
905//________________________________________________________________________
906void AliAnalysisTaskHJetDphi::RunJetQA(const TClonesArray *jetArray, const Double_t rho, THnSparse *hJetPt, THnSparse *hJetArea, THnSparse *hJetQA)
907{
908 Int_t nJets = jetArray->GetEntries();
0330533a 909 Double_t jetPt = 0, jetEta = 0, jetPhi = 0, jetArea = 0, jetPtCorr = 0;
910 Int_t nCons = 0;
aac81e6b 911 for(Int_t ij=0; ij<nJets; ij++)
912 {
913 if(fAnaType==0)
914 {
915 AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(jetArray->At(ij));
916 jetPt = jet->Pt();
917 jetEta = jet->Eta();
918 jetPhi = jet->Phi();
919 jetArea = jet->Area();
920 nCons = jet->GetNumberOfConstituents();
921 }
0330533a 922 else if(fAnaType==1)
aac81e6b 923 {
924 AliAODJet* jet = dynamic_cast<AliAODJet*>(jetArray->At(ij));
925 jetPt = jet->Pt();
926 jetEta = jet->Eta();
927 jetPhi = jet->Phi();
928 jetArea = jet->EffectiveAreaCharged();
929 nCons = jet->GetRefTracks()->GetEntriesFast();
930 }
0330533a 931 else
932 return;
aac81e6b 933 if(!IsGoodJet(jetEta)) continue; // eta cut
934
935 jetPtCorr = jetPt-rho*jetArea;
0330533a 936 Double_t fillPt[] = {jetPtCorr, jetPt, fCentrality, (Double_t)fPtHardBin};
aac81e6b 937 hJetPt->Fill(fillPt);
938
0330533a 939 Double_t fillA[] = {jetPtCorr, jetArea, fCentrality, (Double_t)fPtHardBin};
aac81e6b 940 hJetArea->Fill(fillA);
941
605e0d00 942 Double_t fillQA[] = {jetPtCorr, jetPhi*TMath::RadToDeg(), jetEta, jetArea, (Double_t)nCons, fCentrality, (Double_t)fPtHardBin};
943 hJetQA->Fill(fillQA);
aac81e6b 944 }
945}
946
947
948//________________________________________________________________________
949void AliAnalysisTaskHJetDphi::RunLeadTrkQA()
950{
951 if(fIsEmbedding || fTriggerTrkIndex<0) return;
952 Double_t jetPt = -1, jetPhi = -999;
953
954 if(fAnaType==0)
955 {
956 Int_t nJets = fJetArray->GetEntries();
957 for(Int_t ij=0; ij<nJets; ij++)
958 {
959 AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fJetArray->At(ij));
960 if(!IsGoodJet(jet->Eta())) continue; // eta cut
961 for (Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
962 {
963 if(jet->TrackAt(i)==fTriggerTrkIndex)
964 {
965 jetPt = jet->Pt();
966 jetPhi = jet->Phi();
967 break;
968 }
969 }
970 }
971 }
972 if(fAnaType==1)
973 {
974 jetPt = -1;
975 }
976
977
978 if(jetPt<=0) return;
979 AliVParticle *tt = (AliVParticle*) fTrackArray->At(fTriggerTrkIndex);
980 Double_t fill[] = {tt->Pt(), jetPt, tt->Phi()-jetPhi, tt->Pt()/jetPt, fCentrality};
981 fhLeadTrkQA[fTriggerType]->Fill(fill);
982}
983
984
985//________________________________________________________________________
986void AliAnalysisTaskHJetDphi::StudyKtEffects()
987{
988 if(fAnaType==1) return;
989 if(fTriggerTrkIndex<0) return;
990
991 fKtValue = 0; // dummy
992
0330533a 993
aac81e6b 994 AliPicoTrack *tt = (AliPicoTrack*) fTrackArray->At(fTriggerTrkIndex);
995 Double_t triggerPhi = tt->GetTrackPhiOnEMCal();
996
997 Int_t nJets = fDLJetArray->GetEntries();
aac81e6b 998 for(Int_t ij=0; ij<nJets; ij++)
999 {
1000 AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fDLJetArray->At(ij));
1001 if(!IsGoodJet(jet->Eta())) continue; // eta cut
0330533a 1002 Double_t jetPhi = jet->Phi();
1003 Double_t jetPt = jet->Pt();
1004 Double_t dPhi = CalculateDPhi(triggerPhi,jetPhi);
aac81e6b 1005 if(dPhi<pi+0.6 && dPhi>pi-0.6)
1006 {
0330533a 1007 Double_t fill[] = {tt->Pt(), jetPt, jetPt*TMath::Tan(tt->GetTrackPtOnEMCal()), fCentrality, (Double_t)fPtHardBin};
aac81e6b 1008 fhKtEffects[fTriggerType]->Fill(fill);
1009 }
1010 }
1011}
1012
1013
1014//________________________________________________________________________
1015Bool_t AliAnalysisTaskHJetDphi::RetrieveArraies()
1016{
1017 if(fAnaType==0)
1018 {
1019 // Get mc particles
1020 if (!fMcParticleArrName.IsNull())
1021 {
1022 AliInfo(Form("Retrieve mc particles %s!", fMcParticleArrName.Data()));
1023 fMcParticleArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fMcParticleArrName));
1024 if (!fMcParticleArray)
1025 {
1026 AliError(Form("Could not retrieve mc particles %s!", fMcParticleArrName.Data()));
1027 return kFALSE;
1028 }
1029
1030 TString fMcParticleMapArrName = fMcParticleArrName + "_Map";
1031 fMcParticlelMap = dynamic_cast<AliNamedArrayI*>(fEvent->FindListObject(fMcParticleMapArrName));
1032 if (!fMcParticlelMap)
1033 {
1034 AliWarning(Form("%s: Could not retrieve map for MC particles %s! Will assume MC labels consistent with indexes...", GetName(), fMcParticleArrName.Data()));
1035 return kFALSE;
1036 }
1037 }
1038
1039 // Get track collection
1040 if (!fTrackArrName.IsNull())
1041 {
1042 AliInfo(Form("Retrieve tracks %s!", fTrackArrName.Data()));
1043 fTrackArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTrackArrName));
1044 if (!fTrackArray)
1045 {
1046 AliError(Form("Could not retrieve tracks %s!", fTrackArrName.Data()));
1047 return kFALSE;
1048 }
1049 }
1050
1051 // Get mixed track collection
1052 if (!fEmbTrkArrName.IsNull())
1053 {
1054 AliInfo(Form("Retrieve PYTHIA+PbPb tracks %s!", fEmbTrkArrName.Data()));
1055 fEmbTrkArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fEmbTrkArrName));
1056 if (!fEmbTrkArray)
1057 {
1058 AliError(Form("Could not retrieve PYTHIA+PbPb tracks %s!", fEmbTrkArrName.Data()));
1059 return kFALSE;
1060 }
1061 }
1062
1063 // Get Rho value
1064 if(fCollisionSystem=="PbPb")
1065 {
1066 if(!fRhoName.IsNull())
1067 {
1068 AliInfo(Form("Retrieve rho %s!", fRhoName.Data()));
1069 fRho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoName));
1070 if(!fRho)
1071 {
1072 AliError(Form("Could not retrieve rho %s!",fRhoName.Data()));
1073 return kFALSE;
1074 }
1075 }
1076 }
1077
1078 // Get jet collection
1079 if (!fJetArrName.IsNull())
1080 {
1081 AliInfo(Form("Retrieve jets %s!", fJetArrName.Data()));
1082 fJetArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fJetArrName));
1083 if (!fJetArray)
1084 {
1085 AliError(Form("%s: Could not retrieve jets %s!", GetName(), fJetArrName.Data()));
1086 return kFALSE;
1087 }
1088 }
1089
1090 // Get DL jet collection
1091 if (!fDLJetArrName.IsNull())
1092 {
1093 AliInfo(Form("Retrieve DL jets %s!", fDLJetArrName.Data()));
1094 fDLJetArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fDLJetArrName));
1095 if (!fDLJetArray)
1096 {
1097 AliError(Form("%s: Could not retrieve DL jets %s!", GetName(), fDLJetArrName.Data()));
1098 return kFALSE;
1099 }
1100 }
1101
1102 // Get PL jet collection
1103 if (!fPLJetArrName.IsNull())
1104 {
1105 AliInfo(Form("Retrieve PL jets %s!", fPLJetArrName.Data()));
1106 fPLJetArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fPLJetArrName));
1107 if (!fPLJetArray)
1108 {
1109 AliError(Form("%s: Could not retrieve PL jets %s!", GetName(), fPLJetArrName.Data()));
1110 return kFALSE;
1111 }
1112 }
1113
1114 if(fIsEmbedding)
1115 {
1116 if(fAnaType==0 && !fPtHardBinName)
1117 {
1118 // Get embedded pt hard bin number
1119 fPtHardBinName = static_cast<AliNamedString*>(fEvent->FindListObject("AODEmbeddingFile"));
1120 if(!fPtHardBinName)
1121 {
1122 AliError("The object for pt hard bin information is not available!");
1123 return kFALSE;
1124 }
1125 }
1126 }
1127 }
1128
1129 if(fAnaType==1)
1130 {
1131 // Get mc particles
1132 if (!fMcParticleArrName.IsNull())
1133 {
1134 AliInfo(Form("Retrieve mc particles %s!", fMcParticleArrName.Data()));
1135 if(fAODOut && !fMcParticleArray) fMcParticleArray = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fMcParticleArrName));
1136 if(fAODIn && !fMcParticleArray) fMcParticleArray = dynamic_cast<TClonesArray*>(fAODIn ->FindListObject(fMcParticleArrName));
1137 if(fAODExtension && !fMcParticleArray) fMcParticleArray = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fMcParticleArrName));
1138 if (!fMcParticleArray)
1139 {
1140 AliError(Form("Could not retrieve mc particles %s!", fMcParticleArrName.Data()));
1141 return kFALSE;
1142 }
1143 }
1144
1145 // Get tracks
1146 if (!fTrackArrName.IsNull())
1147 {
1148 AliInfo(Form("Retrieve tracks %s!", fTrackArrName.Data()));
1149 if(fAODIn && !fTrackArray) fTrackArray = dynamic_cast<TClonesArray*>(fAODIn ->FindListObject(fTrackArrName));
1150 if(fAODOut && !fTrackArray) fTrackArray = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fTrackArrName));
1151 if (!fTrackArray)
1152 {
1153 AliError(Form("Could not retrieve tracks %s!", fTrackArrName.Data()));
1154 return kFALSE;
1155 }
1156 }
1157
1158 // Get PYTHIA+PbPb tracks
1159 if (!fEmbTrkArrName.IsNull())
1160 {
1161 AliInfo(Form("Retrieve PYTHIA+PbPb tracks %s!", fEmbTrkArrName.Data()));
1162 if(fAODOut && !fEmbTrkArray) fEmbTrkArray = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fEmbTrkArrName));
1163 if(fAODIn && !fEmbTrkArray) fEmbTrkArray = dynamic_cast<TClonesArray*>(fAODIn ->FindListObject(fEmbTrkArrName));
1164 if(fAODExtension && !fEmbTrkArray) fEmbTrkArray = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fEmbTrkArrName));
1165 if (!fTrackArray)
1166 {
1167 AliError(Form("Could not retrieve PYTHIA+PbPb tracks %s!", fTrackArrName.Data()));
1168 return kFALSE;
1169 }
1170 }
1171
1172 // Get jet collection
1173 if (!fJetArrName.IsNull())
1174 {
1175 AliInfo(Form("Retrieve jets %s!", fJetArrName.Data()));
1176 if(fAODOut && !fJetArray) fJetArray = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetArrName));
1177 if(fAODExtension && !fJetArray) fJetArray = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetArrName));
1178 if(fAODIn && !fJetArray) fJetArray = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetArrName));
1179 if (!fJetArray)
1180 {
1181 AliError(Form("%s: Could not retrieve jets %s!", GetName(), fJetArrName.Data()));
1182 return kFALSE;
1183 }
1184 }
1185
1186 // Get DL jet collection
1187 if (!fDLJetArrName.IsNull())
1188 {
0330533a 1189 AliInfo(Form("Retrieve DL jets %s!", fDLJetArrName.Data()));
aac81e6b 1190 if(fAODOut && !fDLJetArray) fDLJetArray = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fDLJetArrName));
1191 if(fAODExtension && !fDLJetArray) fDLJetArray = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fDLJetArrName));
1192 if(fAODIn && !fDLJetArray) fDLJetArray = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fDLJetArrName));
1193 if (!fDLJetArray)
1194 {
1195 AliError(Form("%s: Could not retrieve DL jets %s!", GetName(), fDLJetArrName.Data()));
1196 return kFALSE;
1197 }
1198 }
1199
1200 // Get PL jet collection
1201 if (!fPLJetArrName.IsNull())
1202 {
1203 AliInfo(Form("Retrieve PL jets %s!", fPLJetArrName.Data()));
1204 if(fAODOut && !fPLJetArray) fPLJetArray = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fPLJetArrName));
1205 if(fAODExtension && !fPLJetArray) fPLJetArray = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fPLJetArrName));
1206 if(fAODIn && !fPLJetArray) fPLJetArray = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fPLJetArrName));
1207 if (!fPLJetArray)
1208 {
1209 AliError(Form("%s: Could not retrieve PL jets %s!", GetName(), fPLJetArrName.Data()));
1210 return kFALSE;
1211 }
1212 }
1213
1214 // Get Rho
1215 if(fCollisionSystem=="PbPb" && !fRhoName.IsNull())
1216 {
1217 AliInfo(Form("Retrieve rho %s!", fRhoName.Data()));
1218 if(fAODOut && !fEvtBkg ) fEvtBkg = (AliAODJetEventBackground*)(fAODOut->FindListObject(fRhoName));
1219 if(fAODExtension && !fEvtBkg ) fEvtBkg = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fRhoName));
1220 if(fAODIn && !fEvtBkg ) fEvtBkg = (AliAODJetEventBackground*)(fAODIn->FindListObject(fRhoName));
1221 if(!fEvtBkg)
1222 {
1223 AliError(Form("Could not retrieve rho %s!",fRhoName.Data()));
1224 return kFALSE;
1225 }
1226 }
1227 }
1228
1229 return kTRUE;
1230}
1231
1232//________________________________________________________________________
1233Double_t AliAnalysisTaskHJetDphi::CalculatePhi(const Double_t py, const Double_t px)
1234{
1235 Double_t phi = TMath::ATan2(py,px);
1236 if(phi<0) phi += 2*pi;
1237 return phi;
1238}
1239
1240//________________________________________________________________________
1241Double_t AliAnalysisTaskHJetDphi::CalculateDPhi(const Double_t phi1, const Double_t phi2)
1242{
1243 Double_t dPhi = phi1-phi2;
1244 if(dPhi>2*pi) dPhi -= 2*pi;
1245 if(dPhi<-2*pi) dPhi += 2*pi;
1246 if(dPhi<-0.5*pi) dPhi += 2*pi;
1247 if(dPhi>1.5*pi) dPhi -= 2*pi;
1248 return dPhi;
1249}
1250
1251//________________________________________________________________________
1252Int_t AliAnalysisTaskHJetDphi::LocateToTPCHole(const Double_t phi)
1253{
1254 if(phi<4-pi/2 && phi>5.5-1.5*pi) return 1; // away-side
1255 else return 0; //near-side
1256}
1257
1258
1259//________________________________________________________________________
1260Bool_t AliAnalysisTaskHJetDphi::AcceptTrack(AliVParticle *track)
1261{
1262 if(track->Pt()<fMinTrkPt || track->Pt()>fMaxTrkPt) return kFALSE;
1263 if(track->Eta()<fMinTrkEta || track->Eta()>fMaxTrkEta) return kFALSE;
1264 if(track->Phi()<fMinTrkPhi || track->Phi()>fMaxTrkPhi) return kFALSE;
1265 return kTRUE;
1266}
1267
1268//________________________________________________________________________
1269Bool_t AliAnalysisTaskHJetDphi::IsGoodAODtrack(AliVParticle *track)
1270{
1271 AliAODTrack *aodtrack = static_cast<AliAODTrack*>(track);
1272 if( fFilterMask>0)
1273 {
1274 if(!aodtrack->TestFilterBit(fFilterMask) ) return kFALSE;
1275 }
1276 else
1277 {
1278 if(!aodtrack->IsHybridGlobalConstrainedGlobal()) return kFALSE;
1279 }
1280 if( fRequireITSRefit && (aodtrack->GetStatus()&AliESDtrack::kITSrefit)==0 ) return kFALSE;
49bf6461 1281 if (fRequireSharedClsCut)
1282 {
1283 Double_t frac = Double_t(aodtrack->GetTPCnclsS())/Double_t(aodtrack->GetTPCncls());
1284 if (frac > 0.4) return kFALSE;
1285 }
aac81e6b 1286 return kTRUE;
1287}
1288
1289//________________________________________________________________________
1290Bool_t AliAnalysisTaskHJetDphi::IsGoodJet(Double_t jetEta)
1291{
1292 Double_t etaCut = (0.9-fRadius>0.5)?0.5:0.9-fRadius;
1293 if(TMath::Abs(jetEta)>etaCut) return kFALSE;
1294 return kTRUE;
1295}
1296
1297//________________________________________________________________________
1298Int_t AliAnalysisTaskHJetDphi::GetPtHardBin(Double_t ptHard)
1299{
1300 const Int_t nBins = 10;
1301 Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
1302 Int_t bin = -1;
1303 while(bin<nBins-1 && binLimits[bin+1]<ptHard)
1304 {
1305 bin++;
1306 }
1307 return bin;
1308}
1309
1310//________________________________________________________________________
1311Int_t AliAnalysisTaskHJetDphi::GetParticleType(Int_t pdg_input)
1312{
1313 Int_t type = -1;
1314 Int_t pdg = TMath::Abs(pdg_input);
1315 if(pdg==211) type = 0;
1316 else if(pdg==321) type = 1;
1317 else if(pdg==2212) type = 2;
1318 else if(pdg==11) type = 3;
1319 else if(pdg>=3122 && pdg<=3334) type = 4;
1320 else type = 9;
1321 return type;
1322}
1323
1324//________________________________________________________________________
1325Double_t AliAnalysisTaskHJetDphi::GetAODTrackPtRes(AliAODTrack *aodtrack)
1326{
1327 Double_t sigma1Pt2 = -1;
1328 Double_t cov[21] = {0,}, pxpypz[3] = {0,}, xyz[3] = {0,};
1329 AliExternalTrackParam *exParam = new AliExternalTrackParam();
1330 aodtrack->GetCovMatrix(cov);
1331 aodtrack->PxPyPz(pxpypz);
1332 aodtrack->GetXYZ(xyz);
1333 exParam->Set(xyz,pxpypz,cov,aodtrack->Charge());
1334 sigma1Pt2 = exParam->GetSigma1Pt2();
1335 delete exParam;
1336 return sigma1Pt2;
1337}
1338
1339//
1340//________________________________________________________________________
1341//
1342void AliAnalysisTaskHJetDphi::PrintConfig()
1343{
1344 const char *decision[2] = {"no","yes"};
1345 printf("\n\n===== h-jet analysis configuration =====\n");
1346 printf("Input event type: %s - %s\n",fCollisionSystem.Data(),fPeriod.Data());
1347 printf("Is this MC data: %s\n",decision[fIsMC]);
1348 printf("Run on particle level: %s\n",decision[fAnalyzeMCTruth]);
1349 printf("Event type selection: %d\n",fOfflineTrgMask);
0330533a 1350 printf("Is embedding? %s\n",decision[fIsEmbedding]);
aac81e6b 1351 printf("Track filter mask: %d\n",fFilterMask);
1352 printf("Require track to have ITS refit? %s\n",decision[fRequireITSRefit]);
49bf6461 1353 printf("Require to cut on fraction of shared TPC clusters? %s\n",decision[fRequireSharedClsCut]);
aac81e6b 1354 printf("Track pt range: %2.2f < pt < %2.2f\n",fMinTrkPt, fMaxTrkPt);
1355 printf("Track eta range: %2.1f < eta < %2.1f\n",fMinTrkEta, fMaxTrkEta);
1356 printf("Track phi range: %2.0f < phi < %2.0f\n",fMinTrkPhi*TMath::RadToDeg(),fMaxTrkPhi*TMath::RadToDeg());
1357 printf("Cut TT away from boundary: %s with distance = %2.2f\n",decision[fCutTPCBoundary],fDistToTPCBoundary);
605e0d00 1358 printf("Avoid TPC holes: %s\n", decision[fSwitchOnAvoidTpcHole]);
aac81e6b 1359 printf("Jet cone size R = %1.1f, and jet pt > %1.0f GeV/c \n",fRadius,fJetPtMin);
1360 printf("Run track QA: %s\n",decision[fRunTrkQA]);
1361 printf("Run jet QA: %s\n",decision[fRunJetQA]);
1362 printf("Run single inclusive h+jet analysis: %s\n",decision[fRunSingleInclHJet]);
1363 printf("TT interval: %2.0f < pt < %2.0f\n",fTTMinPt, fTTMaxPt);
1364 printf("Run leading track QA: %s\n",decision[fRunLeadTrkQA]);
1365 printf("Run kT effects study: %s\n",decision[fStudyKtEffects]);
1366 printf("Run background flow: %s\n",decision[fRunBkgFlow]);
1367 printf("=======================================\n\n");
1368}
1369
1370//________________________________________________________________________
1371Double_t AliAnalysisTaskHJetDphi::GetZ(const Double_t trkPx, const Double_t trkPy, const Double_t trkPz, const Double_t jetPx, const Double_t jetPy, const Double_t jetPz)
1372{
1373 return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/(jetPx*jetPx+jetPy*jetPy+jetPz*jetPz);
1374}
1375
1376//________________________________________________________________________
1377void AliAnalysisTaskHJetDphi::Terminate(Option_t *)
1378{
1379 // Called once at the end of the query
1380}