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