AliPhysics  857879d (857879d)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskHJetDphi.cxx
Go to the documentation of this file.
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"
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"
50 
51 #include <iostream>
52 using std::cout;
53 using std::cerr;
54 using std::endl;
55 
57 
58 const Double_t pi = TMath::Pi();
59 const Double_t kSector = pi/9;
60 
61 //________________________________________________________________________
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),
68  fEsdTrkCut(0x0), fEsdHybCut(0x0), fFilterMask(0), fRequireITSRefit(kTRUE), fRequireSharedClsCut(kTRUE),
69  fIsInit(kFALSE), fNonStdFile(""), fMcParticleArrName(""), fMcParticleArray(0x0), fMcParticlelMap(0x0),
70  fEmbTrkArrName(""), fEmbTrkArray(0x0), fTrackArrName(""), fTrackArray(0x0),
71  fTriggerTrkIndex(-1), fTriggerTrkPt(-1), fSwitchOnAvoidTpcHole(kFALSE), fAvoidTpcHole(0), fCutTPCBoundary(kFALSE), fDistToTPCBoundary(0.),
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 //________________________________________________________________________
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),
127  fEsdTrkCut(0x0), fEsdHybCut(0x0), fFilterMask(0), fRequireITSRefit(kTRUE), fRequireSharedClsCut(kTRUE),
128  fIsInit(kFALSE), fNonStdFile(""), fMcParticleArrName(""), fMcParticleArray(0x0), fMcParticlelMap(0x0),
129  fEmbTrkArrName(""), fEmbTrkArray(0x0), fTrackArrName(""), fTrackArray(0x0),
130  fTriggerTrkIndex(-1), fTriggerTrkPt(-1), fSwitchOnAvoidTpcHole(kFALSE), fAvoidTpcHole(0), fCutTPCBoundary(kFALSE), fDistToTPCBoundary(0.),
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 //________________________________________________________________________
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 //________________________________________________________________________
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;
229  const Int_t nBinsJetqa[dimJetqa] = {nJetPtBins/5, 36, 24, 6, 100, 10, 11};
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 = 4;
235  const Int_t nBinsTT[dimTT] = {nTrkPtBins, 10, 11, 10};
236  const Double_t lowBinTT[dimTT] = {lowTrkPtBin, 0, 0, 0};
237  const Double_t hiBinTT[dimTT] = {upTrkPtBin, 100, 11, 10};
238 
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};
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 
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 
264  fhEventStat = new TH1F("fhEventStat","Event statistics for jet analysis",8,0,8);
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");
272  fhEventStat->GetXaxis()->SetBinLabel(8,"kJetService");
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);
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  {
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);
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 
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 
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);
349  fHJetPhiCorr[i][j]->SetBinEdges(3,newbins);
350  fOutputList->Add(fHJetPhiCorr[i][j]);
351  }
352  }
353 
354 
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 
360 
361  if(fRunBkgFlow)
362  {
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);
364  fOutputList->Add(fHJetPhiCorrUp[i]);
365 
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);
367  fOutputList->Add(fHJetPhiCorrDown[i]);
368  }
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 //________________________________________________________________________
419 {
420  AliInfo("User Nofity");
421 
422  Int_t runNumber = InputEvent()->GetRunNumber();
423 
424  fAvoidTpcHole = 0;
426  {
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 
430  for(Int_t i=0; i<28; i++)
431  {
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])
441  {
442  fAvoidTpcHole = 2;
443  break;
444  }
445  }
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 //________________________________________________________________________
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  {
527  fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
528  }
529  else
530  {
531  if(fAnaType==1) fEvent = AODEvent();
532  }
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 
550  if(fVerbosity>1)
551  {
552  TList *list = 0x0;
553  if(fAnaType==0) list = fEvent->GetList();
554  else list = fAODOut->GetList();
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
563  if(!fIsInit) return;
564 
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;
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  }
612 
614  fhEventStat->Fill(7.5);
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;
626 
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  }
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.8, 0x0, fHJetPhiCorrUp[fTriggerType]);
649  RunSingleInclHJetCorr(trigPt, trigPhi, trigEta, fJetArray, fRhoValue-1.8, 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();
682  }
683 
685 
686  PostData(1, fOutputList);
687  return;
688 }
689 
690 //________________________________________________________________________
691 Int_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 
740  if(trigIndex>-1)
741  {
742  AliVParticle *tt = (AliVParticle*) fTrackArray->At(trigIndex);
743  trigPt = tt->Pt();
744  trigPhi = tt->Phi();
745  trigEta = tt->Eta();
746 
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  }
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  }
762 
763  if(trigIndex==-1) { trigPt = -1; trigPhi = -999; trigEta = -999;}
764  if(arrayType<2) fTriggerTrkIndex = trigIndex;
765  return nTT;
766 }
767 
768 //________________________________________________________________________
769 void 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  {
775  Double_t fillTT[] = {trigPt, fCentrality, (Double_t)fPtHardBin,static_cast<Double_t>(Entry()%10)};
776  hTT->Fill(fillTT);
777  }
778 
779  Int_t nJets = jetArray->GetEntries();
780  Double_t jetPt = 0, jetEta = 0, jetPhi = 0, jetArea = 0;
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  }
791  else if(fAnaType==1)
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  }
799  else
800  return;
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  {
807  Double_t fill[] = {trigPt,jetPtCorr,dPhi,jetArea,fCentrality,trigEta-jetEta, (Double_t)fPtHardBin,static_cast<Double_t>(Entry()%10)};
808  hn->Fill(fill);
809  }
810  }
811 }
812 
813 
814 //________________________________________________________________________
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  {
856  Double_t fillPhiRes[] = {esdtrack->Pt(),part->Phi()-esdtrack->Phi(),(Double_t)type,fCentrality};
857  fhTrkPhiRes[fTriggerType]->Fill(fillPhiRes);
858  }
859  resolution = (part->Pt()-esdtrack->Pt())/part->Pt();
860  }
861  Double_t fillPtRes[] = {esdtrack->Pt(),esdtrack->Pt()*TMath::Sqrt(esdtrack->GetSigma1Pt2()),resolution,(Double_t)type,fCentrality};
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();
893  Double_t fillPhiRes[] = {aodtrack->Pt(),part->Phi()-aodtrack->Phi(),(Double_t)type,fCentrality};
894  fhTrkPhiRes[fTriggerType]->Fill(fillPhiRes);
895  }
896  if(sigma1Pt2>0)
897  {
898  Double_t fillPtRes[5] = {aodtrack->Pt(),aodtrack->Pt()*TMath::Sqrt(sigma1Pt2),resolution,(Double_t)type,fCentrality};
899  fhTrkPtRes[fTriggerType]->Fill(fillPtRes);
900  }
901  }
902  }
903 }
904 
905 //________________________________________________________________________
906 void AliAnalysisTaskHJetDphi::RunJetQA(const TClonesArray *jetArray, const Double_t rho, THnSparse *hJetPt, THnSparse *hJetArea, THnSparse *hJetQA)
907 {
908  Int_t nJets = jetArray->GetEntries();
909  Double_t jetPt = 0, jetEta = 0, jetPhi = 0, jetArea = 0, jetPtCorr = 0;
910  Int_t nCons = 0;
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  }
922  else if(fAnaType==1)
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  }
931  else
932  return;
933  if(!IsGoodJet(jetEta)) continue; // eta cut
934 
935  jetPtCorr = jetPt-rho*jetArea;
936  Double_t fillPt[] = {jetPtCorr, jetPt, fCentrality, (Double_t)fPtHardBin};
937  hJetPt->Fill(fillPt);
938 
939  Double_t fillA[] = {jetPtCorr, jetArea, fCentrality, (Double_t)fPtHardBin};
940  hJetArea->Fill(fillA);
941 
942  Double_t fillQA[] = {jetPtCorr, jetPhi*TMath::RadToDeg(), jetEta, jetArea, (Double_t)nCons, fCentrality, (Double_t)fPtHardBin};
943  hJetQA->Fill(fillQA);
944  }
945 }
946 
947 
948 //________________________________________________________________________
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 //________________________________________________________________________
987 {
988  if(fAnaType==1) return;
989  if(fTriggerTrkIndex<0) return;
990 
991  fKtValue = 0; // dummy
992 
993 
995  Double_t triggerPhi = tt->GetTrackPhiOnEMCal();
996 
997  Int_t nJets = fDLJetArray->GetEntries();
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
1002  Double_t jetPhi = jet->Phi();
1003  Double_t jetPt = jet->Pt();
1004  Double_t dPhi = CalculateDPhi(triggerPhi,jetPhi);
1005  if(dPhi<pi+0.6 && dPhi>pi-0.6)
1006  {
1007  Double_t fill[] = {tt->Pt(), jetPt, jetPt*TMath::Tan(tt->GetTrackPtOnEMCal()), fCentrality, (Double_t)fPtHardBin};
1008  fhKtEffects[fTriggerType]->Fill(fill);
1009  }
1010  }
1011 }
1012 
1013 
1014 //________________________________________________________________________
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  {
1189  AliInfo(Form("Retrieve DL jets %s!", fDLJetArrName.Data()));
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 //________________________________________________________________________
1234 {
1235  Double_t phi = TMath::ATan2(py,px);
1236  if(phi<0) phi += 2*pi;
1237  return phi;
1238 }
1239 
1240 //________________________________________________________________________
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 //________________________________________________________________________
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 //________________________________________________________________________
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 //________________________________________________________________________
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;
1282  {
1283  Double_t frac = Double_t(aodtrack->GetTPCnclsS())/Double_t(aodtrack->GetTPCncls());
1284  if (frac > 0.4) return kFALSE;
1285  }
1286  return kTRUE;
1287 }
1288 
1289 //________________________________________________________________________
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 //________________________________________________________________________
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 //________________________________________________________________________
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 //________________________________________________________________________
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 //
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);
1350  printf("Is embedding? %s\n",decision[fIsEmbedding]);
1351  printf("Track filter mask: %d\n",fFilterMask);
1352  printf("Require track to have ITS refit? %s\n",decision[fRequireITSRefit]);
1353  printf("Require to cut on fraction of shared TPC clusters? %s\n",decision[fRequireSharedClsCut]);
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);
1358  printf("Avoid TPC holes: %s\n", decision[fSwitchOnAvoidTpcHole]);
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 //________________________________________________________________________
1371 Double_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 //________________________________________________________________________
1378 {
1379  // Called once at the end of the query
1380 }
Double_t GetTrackPhiOnEMCal() const
Definition: AliPicoTrack.h:62
Int_t pdg
const Double_t kSector
Double_t Area() const
Definition: AliEmcalJet.h:123
Double_t 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)
double Double_t
Definition: External.C:58
Double_t fCentrality
Trigger type of the event.
Int_t fPtHardBin
Pt hard bin param.
AliESDtrackCuts * fEsdHybCut
Track cuts for global tracks in ESD.
Definition: External.C:236
Double_t Eta() const
Definition: AliEmcalJet.h:114
AliNamedString * fPtHardBinName
Event background for LEGO train.
TString fileName
TClonesArray * fPLJetArray
Array of the found jets.
Double_t Phi() const
Definition: AliEmcalJet.h:110
Bool_t IsGoodAODtrack(AliVParticle *track)
TSystem * gSystem
centrality
void UserExec(Option_t *option)
TList * list
TDirectory file where lists per trigger are stored in train ouput.
Int_t fCollisionSystem
AliAODJetEventBackground * fEvtBkg
Value of the rho parameter.
Bool_t AcceptTrack(AliVParticle *track)
Double_t GetTrackPtOnEMCal() const
Definition: AliPicoTrack.h:64
AliVEvent::EOfflineTriggerTypes fOfflineTrgMask
where we take the jets from can be input or output AOD
Int_t FindSingleIncTrigger(const TClonesArray *trackArray, Double_t &trigPt, Double_t &trigPhi, Double_t &trigEta, const Int_t arrayType)
UInt_t fFilterMask
Track cuts for complementary tracks in ESD.
static Bool_t Selected(Bool_t bSet=kFALSE, Bool_t bNew=kTRUE)
UShort_t GetNumberOfConstituents() const
Definition: AliEmcalJet.h:133
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:153
Bool_t IsGoodJet(Double_t jetEta)
Int_t GetPtHardBin(Double_t ptHard)
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:132
Bool_t fSwitchOnAvoidTpcHole
Trigger track pt.
TString fNonStdFile
Flag if all the arraies are successfully retrieved.
int Int_t
Definition: External.C:63
Int_t GetParticleType(Int_t pdg_input)
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
TString fTrackArrName
Array of PbPb tracks + PYTHIA tracks.
AliAODEvent * fAODIn
ESD event.
Double_t GetAODTrackPtRes(AliAODTrack *track)
void RunJetQA(const TClonesArray *jetArray, const Double_t rho, THnSparse *hJetPt, THnSparse *hJetArea, THnSparse *hJetQA)
Double_t CalculatePhi(const Double_t py, const Double_t px)
Double_t Pt() const
Definition: AliPicoTrack.h:23
TString fRhoName
Array of the embedded PYTHIA jet array on detector level.
AliAODExtension * fAODExtension
Output AOD event.
Double_t Pt() const
Definition: AliEmcalJet.h:102
ClassImp(AliAnalysisTaskHJetDphi) const Double_t pi
TString fEmbTrkArrName
Array of mc map for EMCal train.
Double_t CalculateDPhi(const Double_t phi1, const Double_t phi2)
Int_t fTriggerTrkIndex
Array of input tracks.
AliESDEvent * fESD
Input event.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
TClonesArray * fDLJetArray
Array of the embedded PYTHIA jet array on particle level.
const char Option_t
Definition: External.C:48
AliVEvent * fEvent
MC events.
Int_t fAODfilterBits[2]
Random number generator.
const Double_t pi
bool Bool_t
Definition: External.C:53
AliAODEvent * fAODOut
Input AOD event.
Double_t fRhoValue
Rho parameter.
Int_t LocateToTPCHole(const Double_t phi)
Double_t fMaxVtxZ
V0M for current event.
void RunSingleInclHJetCorr(Double_t trigPt, Double_t trigPhi, Double_t trigEta, const TClonesArray *jetArray, Double_t rho, THnSparse *hTT, THnSparse *hn)
AliNamedArrayI * fMcParticlelMap
Array of input mc particles.
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65