]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskHJetEmbed.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
bug-fix: rotation of sub-leading jet in di-jet
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskHJetEmbed.cxx
CommitLineData
ee87a008 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 <TParameter.h>
18
19#include "AliAODEvent.h"
20#include "AliAODInputHandler.h"
21#include "AliAnalysisManager.h"
22#include "AliAnalysisTask.h"
23#include "AliCentrality.h"
24#include "AliAnalysisTaskHJetEmbed.h"
25#include "AliESDEvent.h"
26#include "AliESDInputHandler.h"
27#include "AliVParticle.h"
28#include "AliVTrack.h"
29#include "AliInputEventHandler.h"
30#include "AliMCEvent.h"
31#include "AliStack.h"
32#include "AliGenEventHeader.h"
33#include "AliGenPythiaEventHeader.h"
34#include "AliLog.h"
35#include "AliRhoParameter.h"
aaea4d43 36#include "AliNamedString.h"
ee87a008 37
38#include "AliEmcalJet.h"
39
40#include <iostream>
41using std::cout;
42using std::cerr;
43using std::endl;
44
45ClassImp(AliAnalysisTaskHJetEmbed)
46
47const Double_t pi = TMath::Pi();
d81a6f29 48//const Double_t areaCut[4] = {0.1, 0.23, 0.4, 0.63};
ee87a008 49
50//________________________________________________________________________
51AliAnalysisTaskHJetEmbed::AliAnalysisTaskHJetEmbed() :
52 AliAnalysisTaskSE(),
53 fVerbosity(0), fAnaType(1), fPeriod("lhc11h"), fCollisionSystem("PbPb"),
54 fEvent(0), fTriggerType(-1), fCentrality(-1), fMaxVtxZ(10),
55 fMCParticleArrName(""), fMCParticleArray(0x0),
56 fTrackArrName(""), fTrackArray(0x0), fTriggerTrkIndex(-1),
57 fMinTrkPt(0.15), fMaxTrkPt(1e4), fMinTrkEta(-0.9), fMaxTrkEta(0.9), fMinTrkPhi(0), fMaxTrkPhi(2*pi),
aaea4d43 58 fTTtype(0),
ee87a008 59 fRadius(0.4), fJetArrName(""), fPLJetArrName(""), fDLJetArrName(""),
60 fJetArray(0x0), fPLJetArray(0x0), fDLJetArray(0x0),
61 fRhoName(""), fRho(0x0), fRhoValue(0),
aaea4d43 62 fPtHardBinName(0x0), fPtHardBin(-1), fRandom(0x0),
ee87a008 63 fRunQA(kTRUE), fRunHJet(kTRUE), fRunMatch(kTRUE),
aaea4d43 64 fRunPL(kFALSE), fRunDL(kTRUE),
ee87a008 65 fOutputList(0x0), fhEventStat(0x0), fhPtHardBins(0x0)
66{
67 // Constructor
68
69 // Define input and output slots here
70 // Input slot #0 works with a TChain
71 //DefineInput(0, TChain::Class());
72 //DefineOutput(1, TList::Class());
73 // Output slot #0 id reserved by the base class for AOD
74
75 for(Int_t i=0; i<kNTrig; i++)
76 {
77 fhVtxZ[i] = 0x0;
78 fhCentrality[i] = 0x0;
79 fhRhoVsCent[i] = 0x0;
80
81 fhDLJetPtVsCent[i] = 0x0;
82 fhPLJetPtVsCent[i] = 0x0;
83
84 fhPLTT[i] = 0x0;
85 fhDLTT[i] = 0x0;
86 fhPLHJet[i] = 0x0;
87 fhDLHJet[i] = 0x0;
88 fhTTPtQA[i] = 0x0;
89 fhTTPt[i] = 0x0;
90 fhHJet[i] = 0x0;
91
92 fhJetPtGeoMatch[i] = 0x0;
93 fhJetPtEnMatch[i] = 0x0;
94 fhJetPhiGeoMatch[i] = 0x0;
95 fhJetPhiEnMatch[i] = 0x0;
96 }
aaea4d43 97 for(Int_t i=0; i<kNTT; i++)
98 {
99 if(i==0)
100 { fMinTTPt[i] = 19; fMaxTTPt[i] = 25; }
101 else
102 { fMinTTPt[i] = -1; fMaxTTPt[i] = -1; }
103 }
ee87a008 104}
105
106//________________________________________________________________________
107AliAnalysisTaskHJetEmbed::AliAnalysisTaskHJetEmbed(const char *name) :
108 AliAnalysisTaskSE(name),
109 fVerbosity(0), fAnaType(1), fPeriod("lhc11h"), fCollisionSystem("PbPb"),
110 fEvent(0), fTriggerType(-1), fCentrality(-1), fMaxVtxZ(10),
111 fMCParticleArrName(""), fMCParticleArray(0x0),
112 fTrackArrName(""), fTrackArray(0x0), fTriggerTrkIndex(-1),
113 fMinTrkPt(0.15), fMaxTrkPt(1e4), fMinTrkEta(-0.9), fMaxTrkEta(0.9), fMinTrkPhi(0), fMaxTrkPhi(2*pi),
aaea4d43 114 fTTtype(0),
ee87a008 115 fRadius(0.4), fJetArrName(""), fPLJetArrName(""), fDLJetArrName(""),
116 fJetArray(0x0), fPLJetArray(0x0), fDLJetArray(0x0),
117 fRhoName(""), fRho(0x0), fRhoValue(0),
aaea4d43 118 fPtHardBinName(0x0), fPtHardBin(-1), fRandom(0x0),
ee87a008 119 fRunQA(kTRUE), fRunHJet(kTRUE), fRunMatch(kTRUE),
aaea4d43 120 fRunPL(kFALSE), fRunDL(kTRUE),
ee87a008 121 fOutputList(0x0), fhEventStat(0x0), fhPtHardBins(0x0)
122{
123 // Constructor
124
125 DefineInput(0, TChain::Class());
126 DefineOutput(1, TList::Class());
127
128 for(Int_t i=0; i<kNTrig; i++)
129 {
130 fhVtxZ[i] = 0x0;
131 fhCentrality[i] = 0x0;
132 fhRhoVsCent[i] = 0x0;
133
134 fhDLJetPtVsCent[i] = 0x0;
135 fhPLJetPtVsCent[i] = 0x0;
136
137 fhPLTT[i] = 0x0;
138 fhDLTT[i] = 0x0;
139 fhPLHJet[i] = 0x0;
140 fhDLHJet[i] = 0x0;
141 fhTTPtQA[i] = 0x0;
142 fhTTPt[i] = 0x0;
143 fhHJet[i] = 0x0;
144
145 fhJetPtGeoMatch[i] = 0x0;
146 fhJetPtEnMatch[i] = 0x0;
147 fhJetPhiGeoMatch[i] = 0x0;
148 fhJetPhiEnMatch[i] = 0x0;
149 }
aaea4d43 150 for(Int_t i=0; i<kNTT; i++)
151 {
152 if(i==0)
153 { fMinTTPt[i] = 19; fMaxTTPt[i] = 25; }
154 else
155 { fMinTTPt[i] = -1; fMaxTTPt[i] = -1; }
156 }
ee87a008 157}
158//________________________________________________________________________
159AliAnalysisTaskHJetEmbed::~AliAnalysisTaskHJetEmbed()
160{
161 //Destructor
162 if(fRho) delete fRho;
163 if(fOutputList) { fOutputList->Delete(); delete fOutputList;}
164}
165
166//________________________________________________________________________
167void AliAnalysisTaskHJetEmbed::UserCreateOutputObjects()
168{
169 // Create histograms
170
171 const Int_t nJetPtBins = 40;
172 const Float_t lowJetPtBin=-50, upJetPtBin=150;
173
174 const Int_t nTrkPtBins = 100;
175 const Float_t lowTrkPtBin=0, upTrkPtBin=100;
176
177
178 // QA
179 const Int_t dimJetqa = 3;
180 const Int_t nBinsJetqa[dimJetqa] = {nJetPtBins, 30, 11};
181 const Double_t lowBinJetqa[dimJetqa] = {lowJetPtBin, 0, 0};
182 const Double_t hiBinJetqa[dimJetqa] = {upJetPtBin, 30, 11};
183
184 // h+jet
ee87a008 185 const Int_t dimTT = 3;
186 const Int_t nBinsTT[dimTT] = {nTrkPtBins, 30, 11};
187 const Double_t lowBinTT[dimTT] = {lowTrkPtBin, 0, 0};
188 const Double_t hiBinTT[dimTT] = {upTrkPtBin, 30, 11};
189
190 const Int_t dimHJet = 6;
3ec095a2 191 const Int_t nBinsHJet[dimHJet] = {nTrkPtBins, nJetPtBins, 140, 8, 30, 11};
192 const Double_t lowBinHJet[dimHJet] = {lowTrkPtBin, lowJetPtBin, pi-4.95, 0, 0, 0};
193 const Double_t hiBinHJet[dimHJet] = {upTrkPtBin, upJetPtBin, pi+2.05, 0.8, 30, 11};
ee87a008 194
195 // Match
196 const Int_t dimMthPt = 6;
197 const Int_t nBinsMthPt[dimMthPt] = {nJetPtBins, nJetPtBins, 20, 20, 30, 11};
198 const Double_t lowBinMthPt[dimMthPt] = {lowJetPtBin, lowJetPtBin, -0.95, 0, 0, 0};
199 const Double_t hiBinMthPt[dimMthPt] = {upJetPtBin, upJetPtBin, 1.05, 1, 30, 11};
200
aaea4d43 201 const Int_t dimMthPhi = 7;
202 const Int_t nBinsMthPhi[dimMthPhi] = {nTrkPtBins, nJetPtBins/2, 70, 70, 10, 1, 11};
203 const Double_t lowBinMthPhi[dimMthPhi] = {lowTrkPtBin, lowJetPtBin, pi-4.95, pi-4.95, 0, 0, 0};
204 const Double_t hiBinMthPhi[dimMthPhi] = {upTrkPtBin, upJetPtBin, pi+2.05, pi+2.05, 1, 10, 11};
ee87a008 205
206 OpenFile(1);
207 fOutputList = new TList();
208 fOutputList->SetOwner(1);
209
210 fhEventStat = new TH1F("fhEventStat","Event statistics for jet analysis",9,0,9);
211 fhEventStat->GetXaxis()->SetBinLabel(1,"All");
212 fhEventStat->GetXaxis()->SetBinLabel(2,"PS");
213 fhEventStat->GetXaxis()->SetBinLabel(3,"Vtx");
214 fhEventStat->GetXaxis()->SetBinLabel(4,"Vtx+10cm");
215 fhEventStat->GetXaxis()->SetBinLabel(5,"kMB");
216 fhEventStat->GetXaxis()->SetBinLabel(6,"kCentral");
217 fhEventStat->GetXaxis()->SetBinLabel(7,"kSemiCentral");
218 fhEventStat->GetXaxis()->SetBinLabel(8,"kEMCEGA");
219 fhEventStat->GetXaxis()->SetBinLabel(9,"kEMCEJE");
220 fOutputList->Add(fhEventStat);
221
222 fhPtHardBins = new TH1F("fhPtHardBins","Number of events in each pT hard bin",11,0,11);
223 fOutputList->Add(fhPtHardBins);
224
225 const char *triggerName[kNTrig] = {"kMB","kEGA","kEJE"};
226
227 for(Int_t i=0; i<kNTrig; i++)
228 {
229 if( i!=0 ) continue;
230
231 fhVtxZ[i] = new TH1F(Form("%s_fhVtxZ",triggerName[i]),Form("%s: z distribution of event vertexz;z(cm)",triggerName[i]),400,-20,20);
232 fOutputList->Add(fhVtxZ[i]);
233
234 fhCentrality[i] = new TH1F(Form("%s_fhCentrality",triggerName[i]),Form("%s: Event centrality;centrality",triggerName[i]),100,0,100);
235 fOutputList->Add(fhCentrality[i]);
236
237 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);
238 fOutputList->Add(fhRhoVsCent[i]);
239
240 if(fRunQA)
241 {
242 fhPLJetPtVsCent[i] = new THnSparseF(Form("%s_fhPLJetPtVsCent",triggerName[i]),Form("PYTHIA: jet p_{T} vs centrality vs pT hard bin (particle-level,R=%1.1f);p_{T,jet}^{ch} (GeV/c);centrality;pT hard bin",fRadius),dimJetqa,nBinsJetqa,lowBinJetqa,hiBinJetqa);
243 fOutputList->Add(fhPLJetPtVsCent[i]);
244
245 fhDLJetPtVsCent[i] = new THnSparseF(Form("%s_fhDLJetPtVsCent",triggerName[i]),Form("PYTHIA: jet p_{T} vs centrality vs pT hard bin (detector-level,R=%1.1f);p_{T,jet}^{ch} (GeV/c);centrality;pT hard bin",fRadius),dimJetqa,nBinsJetqa,lowBinJetqa,hiBinJetqa);
246 fOutputList->Add(fhDLJetPtVsCent[i]);
247 }
248
249 if(fRunHJet)
250 {
251 fhPLTT[i] = new THnSparseF(Form("%s_fhPLTT",triggerName[i]),Form("PYTHIA: TT p_{T} vs centrality vs pT hard bin (particle-level);p_{T,TT}^{ch} (GeV/c);centrality;pT hard bin"),dimTT,nBinsTT,lowBinTT,hiBinTT);
252 fOutputList->Add(fhPLTT[i]);
253
254 fhDLTT[i] = new THnSparseF(Form("%s_fhDLTT",triggerName[i]),Form("PYTHIA: TT p_{T} vs centrality vs pT hard bin (detector-level);p_{T,TT}^{ch} (GeV/c);centrality;pT hard bin"),dimTT,nBinsTT,lowBinTT,hiBinTT);
255 fOutputList->Add(fhDLTT[i]);
256
257 fhTTPt[i] = new THnSparseF(Form("%s_fhTTPt",triggerName[i]),Form("Embedded: TT p_{T} vs centrality vs pT hard bin;p_{T,TT}^{ch} (GeV/c);centrality;pT hard bin"),dimTT,nBinsTT,lowBinTT,hiBinTT);
258 fOutputList->Add(fhTTPt[i]);
259
ee87a008 260 fhPLHJet[i] = new THnSparseF(Form("%s_fhPLHJet",triggerName[i]),Form("PYTHIA: TT p_{T} vs jet p_{T} vs #Delta#varphi vs jet area vs centrality vs pT hard bin (particle-level, R=%1.1f);p_{T,TT}^{ch} (GeV/c);p_{T,jet}^{ch} (GeV/c);#Delta#varphi;Area;centrality;pT hard bin",fRadius),dimHJet,nBinsHJet,lowBinHJet,hiBinHJet);
261 fOutputList->Add(fhPLHJet[i]);
262
263 fhDLHJet[i] = new THnSparseF(Form("%s_fhDLHJet",triggerName[i]),Form("PYTHIA: TT p_{T} vs jet p_{T} vs #Delta#varphi vs jet area vs centrality vs pT hard bin (detector-level, R=%1.1f);p_{T,TT}^{ch} (GeV/c);p_{T,jet}^{ch} (GeV/c);#Delta#varphi;Area;centrality;pT hard bin",fRadius),dimHJet,nBinsHJet,lowBinHJet,hiBinHJet);
264 fOutputList->Add(fhDLHJet[i]);
265
266 fhHJet[i] = new THnSparseF(Form("%s_fhHJet",triggerName[i]),Form("Embedded: TT p_{T} vs jet p_{T} vs #Delta#varphi vs jet area vs centrality vs pT hard bin (R=%1.1f);p_{T,TT}^{ch} (GeV/c);p_{T,jet}^{ch} (GeV/c);#Delta#varphi;Area;centrality;pT hard bin",fRadius),dimHJet,nBinsHJet,lowBinHJet,hiBinHJet);
267 fOutputList->Add(fhHJet[i]);
268 }
269
270 if(fRunMatch)
271 {
272 fhJetPtGeoMatch[i] = new THnSparseF(Form("%s_fhJetPtGeoMatch",triggerName[i]),Form("Embed: generated p_{T,jet} vs reconstructed p_{T,jet} vs jet p_{T} difference vs dR vs centrality vs pT hard bin (R=%1.1f);p_{T,jet}^{gen} (GeV/c);p_{T,jet}^{rec} (GeV/c);(p_{T,jet}^{rec}-p_{T,jet}^{gen})/p_{T,jet}^{gen};dR;centrality;pT hard bin",fRadius),dimMthPt,nBinsMthPt,lowBinMthPt,hiBinMthPt);
273 fOutputList->Add(fhJetPtGeoMatch[i]);
274
275 fhJetPtEnMatch[i] = new THnSparseF(Form("%s_fhJetPtEnMatch",triggerName[i]),Form("Embed: generated p_{T,jet} vs reconstructed p_{T,jet} vs jet p_{T} difference vs dR vs centrality vs pT hard bin (R=%1.1f);p_{T,jet}^{gen} (GeV/c);p_{T,jet}^{rec} (GeV/c);(p_{T,jet}^{rec}-p_{T,jet}^{gen})/p_{T,jet}^{gen};dR;centrality;pT hard bin",fRadius),dimMthPt,nBinsMthPt,lowBinMthPt,hiBinMthPt);
276 fOutputList->Add(fhJetPtEnMatch[i]);
277
aaea4d43 278 fhJetPhiGeoMatch[i] = new THnSparseF(Form("%s_fhJetPhiGeoMatch",triggerName[i]),Form("Embed: p_{T,TT} vs p_{T,jet}^{det} vs #Delta#varphi_{TT} vs #Delta#varphi_{jet} vs dR vs centrality vs pT hard bin (R=%1.1f);p_{T,TT} (GeV/c);p_{T,jet}^{det} (GeV/c);#Delta#varphi_{TT};#Delta#varphi_{jet};dR;centrality;pThard bin",fRadius),dimMthPhi,nBinsMthPhi,lowBinMthPhi,hiBinMthPhi);
ee87a008 279 fOutputList->Add(fhJetPhiGeoMatch[i]);
280
aaea4d43 281 fhJetPhiEnMatch[i] = new THnSparseF(Form("%s_fhJetPhiEnMatch",triggerName[i]),Form("Embed: p_{T,TT} vs p_{T,jet}^{det} vs #Delta#varphi_{TT} vs #Delta#varphi_{jet} vs dR vs centrality vs pT hard bin (R=%1.1f);p_{T,TT} (GeV/c);p_{T,jet}^{det} (GeV/c);#Delta#varphi_{TT};#Delta#varphi_{jet};dR;centrality;pThard bin",fRadius),dimMthPhi,nBinsMthPhi,lowBinMthPhi,hiBinMthPhi);
ee87a008 282 fOutputList->Add(fhJetPhiEnMatch[i]);
283 }
284 }
285
286 //error calculation in THnSparse
287 Int_t nObj = fOutputList->GetEntries();
288 for(Int_t i=0; i<nObj; i++)
289 {
290 TObject *obj = (TObject*) fOutputList->At(i);
291 if (obj->IsA()->InheritsFrom( "THnSparse" ))
292 {
293 THnSparseF *hn = (THnSparseF*)obj;
294 hn->Sumw2();
295 }
296 }
297
3ec095a2 298 fRandom = new TRandom3();
299
ee87a008 300 PrintConfig();
301 PostData(1, fOutputList);
302}
303
304//________________________________________________________________________
305void AliAnalysisTaskHJetEmbed::UserExec(Option_t *)
306{
307 // Main loop, called for each event.
308
aaea4d43 309 AliDebug(5,"Entering UserExec");
ee87a008 310 fTriggerType = -1;
ee87a008 311 fEvent = InputEvent();
312 if (!fEvent)
313 {
314 AliError("Input event not available");
315 return;
316 }
aaea4d43 317 AliDebug(5,"Got the input event");
318 if(!fPtHardBinName)
ee87a008 319 {
320 // Get embedded pt hard bin number
aaea4d43 321 fPtHardBinName = static_cast<AliNamedString*>(fEvent->FindListObject("AODEmbeddingFile"));
322 if(!fPtHardBinName)
323 {
324 AliError("The object for pt hard bin information is not available!");
325 return;
326 }
327 }
328 TString fileName = fPtHardBinName->GetString();
329 fileName.Remove(0,50);
330 fileName.Remove(fileName.Index("/"));
331 fPtHardBin = fileName.Atoi();
ee87a008 332 fhEventStat->Fill(0.5);
333 UInt_t trigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
aaea4d43 334
ee87a008 335 if(fPeriod.Contains("lhc11h",TString::kIgnoreCase))
336 {
337 if (trigger & AliVEvent::kAnyINT) { fTriggerType=0; }
338 else if (trigger & AliVEvent::kCentral) { fTriggerType=0; }
339 else if (trigger & AliVEvent::kSemiCentral) { fTriggerType=0; }
340 else if (trigger & AliVEvent::kEMCEGA) { fTriggerType=1; }
341 else if (trigger & AliVEvent::kEMCEJE) { fTriggerType=2; }
342 }
343 else if(fPeriod.Contains("lhc10h",TString::kIgnoreCase))
344 {
345 if (trigger & AliVEvent::kAnyINT) { fTriggerType=0; }
346 }
347 else if(fPeriod.Contains("lhc12a15a",TString::kIgnoreCase))
348 {
349 fTriggerType=0;
350 }
351
352 if(fTriggerType==-1) return;
353 fhEventStat->Fill(1.5);
354
355 // Vertex cut
356 const AliVVertex* vtx = fEvent->GetPrimaryVertex();
357 if (!vtx || vtx->GetNContributors()<1) return;
358 fhEventStat->Fill(2.5);
359 fhVtxZ[fTriggerType]->Fill(vtx->GetZ());
360 if (TMath::Abs(vtx->GetZ())>fMaxVtxZ) return;
361 fhEventStat->Fill(3.5);
362
363 if(fTriggerType==0)
364 {
365 if(trigger & AliVEvent::kCentral) fhEventStat->Fill(5.5);
366 else if (trigger & AliVEvent::kCentral) fhEventStat->Fill(6.5);
367 else fhEventStat->Fill(4.5);
368 }
369 if(fTriggerType==1) fhEventStat->Fill(7.5);
370 if(fTriggerType==2) fhEventStat->Fill(8.5);
371
372 // GetCentrality
373 if(fCollisionSystem=="PbPb")
374 {
375 AliCentrality *centrality = fEvent->GetCentrality();
376 if (centrality)
377 fCentrality = centrality->GetCentralityPercentile("V0M");
378 else
379 fCentrality = 99;
380 }
381 else if(fCollisionSystem=="pp")
382 {
383 fCentrality = 0;
384 }
385
386 // Get track collection and run QA
387 if (!fTrackArray)
388 {
389 fTrackArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTrackArrName));
390 if (!fTrackArray)
391 {
392 AliError(Form("Could not retrieve tracks %s!", fTrackArrName.Data()));
393 return;
394 }
395 if (!fTrackArray->GetClass()->GetBaseClass("AliVParticle"))
396 {
397 AliError(Form("%s: Collection %s does not contain AliVParticle objects!", GetName(), fTrackArrName.Data()));
398 fTrackArray = 0;
399 return;
400 }
401 }
402
ee87a008 403 // Get MC particle array
aaea4d43 404 if (fRunPL && !fMCParticleArray)
ee87a008 405 {
406 fMCParticleArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fMCParticleArrName));
407 if (!fMCParticleArray)
408 {
409 AliError(Form("Could not retrieve tracks %s!", fMCParticleArrName.Data()));
410 return;
411 }
412 }
413
414 // Get Rho value
415 if(fCollisionSystem=="PbPb")
416 {
417 if(!fRho && !fRhoName.IsNull())
418 {
419 fRho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoName));
420 if(!fRho)
421 {
422 AliError(Form("Could not retrieve rho %s!",fRhoName.Data()));
423 return;
424 }
425 }
426 if(fRho) fRhoValue = fRho->GetVal();
427 }
428 else if(fCollisionSystem=="pp")
429 {
430 fRhoValue = 0;
431 }
432
433 // Get jet collection
434 if (!fJetArray && !fJetArrName.IsNull())
435 {
436 fJetArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fJetArrName));
437 if (!fJetArray)
438 {
439 AliError(Form("%s: Could not retrieve jets %s!", GetName(), fJetArrName.Data()));
440 return;
441 }
442 if (!fJetArray->GetClass()->GetBaseClass("AliEmcalJet"))
443 {
444 AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJetArrName.Data()));
445 fJetArray = 0;
446 return;
447 }
448 }
449 // Get particle-level jet array
aaea4d43 450 if (fRunPL && !fPLJetArray && !fPLJetArrName.IsNull())
ee87a008 451 {
452 fPLJetArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fPLJetArrName));
453 if (!fPLJetArray)
454 {
455 AliError(Form("%s: Could not retrieve jets %s!", GetName(), fPLJetArrName.Data()));
456 return;
457 }
458 if (!fPLJetArray->GetClass()->GetBaseClass("AliEmcalJet"))
459 {
460 AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fPLJetArrName.Data()));
461 fPLJetArray = 0;
462 return;
463 }
464 }
465 // Get detector-level jet array
aaea4d43 466 if (fRunDL && !fDLJetArray && !fDLJetArrName.IsNull())
ee87a008 467 {
468 fDLJetArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fDLJetArrName));
469 if (!fDLJetArray)
470 {
471 AliError(Form("%s: Could not retrieve jets %s!", GetName(), fDLJetArrName.Data()));
472 return;
473 }
474 if (!fDLJetArray->GetClass()->GetBaseClass("AliEmcalJet"))
475 {
476 AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fDLJetArrName.Data()));
477 fDLJetArray = 0;
478 return;
479 }
480 }
481
482 fhCentrality[fTriggerType]->Fill(fCentrality);
483 fhRhoVsCent[fTriggerType]->Fill(fCentrality,fRhoValue);
484 fhPtHardBins->Fill(fPtHardBin);
485
486 if(fRunQA) RunQA();
aaea4d43 487 if(fRunHJet)
488 {
489 for(Int_t i=0; i<kNTT; i++)
490 RunHJet(fMinTTPt[i],fMaxTTPt[i]);
491 }
ee87a008 492
493 PostData(1, fOutputList);
494 return;
495}
496
497
498//________________________________________________________________________
aaea4d43 499void AliAnalysisTaskHJetEmbed::RunHJet(const Double_t minPt, const Double_t maxPt)
ee87a008 500{
3ec095a2 501 TArrayI arr;
3ec095a2 502 Int_t counter = 0;
aaea4d43 503 Int_t indexPL = -1;
504
505 if(fRunPL)
ee87a008 506 {
aaea4d43 507 // Find trigger track on particle level
508 const Int_t nParticles = fMCParticleArray->GetEntries();
509 Double_t maxPLPt = -1;
510 arr.Set(nParticles);
511 counter = 0;
512 for(Int_t iPart=0; iPart<nParticles; iPart++)
3ec095a2 513 {
aaea4d43 514 AliVParticle *t = static_cast<AliVParticle*>(fMCParticleArray->At(iPart));
515 //if(!t || t->Charge()==0) continue;
516 //if(!AcceptTrack(t)) continue;
517 Double_t pt = t->Pt();
518 if(fTTtype==0) // single inclusive triggers
519 {
520 if (pt<maxPt && pt>=minPt)
521 {
522 arr.AddAt(iPart,counter);
523 counter++;
524 }
525 }
526 else if(fTTtype==1) // leading triggers
3ec095a2 527 {
aaea4d43 528 if(maxPLPt<pt)
529 {
530 maxPLPt = pt;
531 indexPL = iPart;
532 }
3ec095a2 533 }
534 }
aaea4d43 535 arr.Set(counter);
536 if(fTTtype==0)
ee87a008 537 {
aaea4d43 538 if(counter==0) indexPL = -1;
539 else if(counter==1) indexPL = arr.At(0);
540 else
3ec095a2 541 {
aaea4d43 542 Double_t pro = fRandom->Uniform() * counter;
543 indexPL = arr.At(TMath::FloorNint(pro));
3ec095a2 544 }
ee87a008 545 }
aaea4d43 546 arr.Reset();
ee87a008 547 }
3ec095a2 548
ee87a008 549
550 // Find trigger track on detector level and after embedding
551 const Int_t Ntracks = fTrackArray->GetEntries();
3ec095a2 552 Double_t maxDLPt = 0;
553 Int_t indexDL = -1;
554 arr.Set(Ntracks);
555 counter = 0;
ee87a008 556 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks)
557 {
558 AliVParticle *t = static_cast<AliVParticle*>(fTrackArray->At(iTracks));
559 if(!t || t->Charge()==0) continue;
560 if(!AcceptTrack(t)) continue;
ee87a008 561 if(t->GetLabel()!=0)
562 {
aaea4d43 563 //cout<<iTracks<<" "<<t->Pt()<<" "<<t->GetLabel()<<endl;
3ec095a2 564 Double_t pt = t->Pt();
565 if(fTTtype==0)
566 {
aaea4d43 567 if (pt<maxPt && pt>=minPt)
3ec095a2 568 {
569 arr.AddAt(iTracks,counter);
570 counter++;
571 }
572 }
573 else if(fTTtype==1)
ee87a008 574 {
3ec095a2 575 if(maxDLPt<pt)
576 {
577 maxDLPt = pt;
578 indexDL = iTracks;
579 }
ee87a008 580 }
581 }
3ec095a2 582 }
583 arr.Set(counter);
584 if(fTTtype==0)
585 {
586 if(counter==0) indexDL = -1;
587 else if(counter==1) indexDL = arr.At(0);
ee87a008 588 else
589 {
3ec095a2 590 Double_t pro = fRandom->Uniform() * counter;
591 indexDL = arr.At(TMath::FloorNint(pro));
ee87a008 592 }
593 }
3ec095a2 594 arr.Reset();
595
aaea4d43 596 AliDebug(2,Form("TT indices: PL=%d, DL=%d\n",indexPL,indexDL));
ee87a008 597
598 // Run h+jet
aaea4d43 599 if(fRunPL) FillHJetCor(fMCParticleArray, indexPL, fPLJetArray, fhPLTT[fTriggerType], fhPLHJet[fTriggerType], kFALSE);
600 if(fRunDL) FillHJetCor(fTrackArray, indexDL, fDLJetArray, fhDLTT[fTriggerType], fhDLHJet[fTriggerType], kFALSE);
3ec095a2 601 FillHJetCor(fTrackArray, indexDL, fJetArray, fhTTPt[fTriggerType], fhHJet[fTriggerType], kTRUE);
aaea4d43 602
603 if(fRunMatch) RunMatch(fTrackArray, indexDL);
ee87a008 604}
605
606//________________________________________________________________________
3ec095a2 607void AliAnalysisTaskHJetEmbed::FillHJetCor(const TClonesArray *tracks, const Int_t leadingIndex, const TClonesArray *jetArray, THnSparse *hTT, THnSparse *hn, Bool_t isBkg)
ee87a008 608{
609 if(leadingIndex<0) return;
610
611 AliVParticle *tt = (AliVParticle*) tracks->At(leadingIndex);
612 Double_t triggerPt = tt->Pt();
2942f542 613 Double_t fill1[] = {triggerPt, fCentrality, static_cast<Double_t>(fPtHardBin) };
3ec095a2 614 hTT->Fill(fill1);
aaea4d43 615 AliDebug(2,Form("Found a trigger with pt = %2.2f",triggerPt));
3ec095a2 616
ee87a008 617 Double_t triggerPhi = tt->Phi();
618 if(triggerPhi<0) triggerPhi += 2*pi;
ee87a008 619 Int_t nJets = jetArray->GetEntries();
620 for(Int_t ij=0; ij<nJets; ij++)
621 {
622 AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(jetArray->At(ij));
aaea4d43 623 if(!jet) continue;
624 if(!IsGoodJet(jet)) continue; // eta cut
ee87a008 625 Double_t jetPhi = jet->Phi();
626 Double_t jetPt = jet->Pt();
627 Double_t jetArea = jet->Area();
628 Double_t dPhi = CalculateDPhi(triggerPhi,jetPhi);
2942f542 629 Double_t fill[] = {triggerPt,jetPt-jetArea*fRhoValue,dPhi,jetArea,fCentrality,static_cast<Double_t>(fPtHardBin)};
ee87a008 630 if(!isBkg) fill[1] = jetPt;
aaea4d43 631 AliDebug(10,"Fill the histograms");
ee87a008 632 hn->Fill(fill);
633 }
634}
635
636//________________________________________________________________________
aaea4d43 637void AliAnalysisTaskHJetEmbed::RunMatch(const TClonesArray *tracks, const Int_t leadingIndex)
ee87a008 638{
aaea4d43 639 if(leadingIndex<0) return;
640
641 if(!fDLJetArray || !fJetArray)
ee87a008 642 {
643 AliWarning("Jet array is not available.");
644 return;
645 }
aaea4d43 646
647 AliVParticle *tt = (AliVParticle*) tracks->At(leadingIndex);
648 Double_t dR = 999, fraction = -1;
649 Int_t nJets = fDLJetArray->GetEntries();
ee87a008 650 for(Int_t ij=0; ij<nJets; ij++)
651 {
aaea4d43 652 AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fDLJetArray->At(ij));
653 if(!jet) continue;
ee87a008 654 if(!IsGoodJet(jet)) continue; // eta cut
655 Double_t jetPt = jet->Pt();
aaea4d43 656 if(jetPt<10) continue;
ee87a008 657
aaea4d43 658 // energy matching
659 Int_t mthJetIndexEn = FindEnergyMatchedJet(jet,fJetArray,dR,fraction);
660 if(mthJetIndexEn>-1 && fraction>0.5)
ee87a008 661 {
aaea4d43 662 AliEmcalJet* jetMthEn = dynamic_cast<AliEmcalJet*>(fJetArray->At(mthJetIndexEn));
663 if(jetMthEn)
ee87a008 664 {
2942f542 665 Double_t fill[] = {tt->Pt(),jetPt,CalculateDPhi(tt->Phi(),jet->Phi()),CalculateDPhi(jetMthEn->Phi(),jet->Phi()),dR,fCentrality,static_cast<Double_t>(fPtHardBin)};
aaea4d43 666 fhJetPhiEnMatch[fTriggerType]->Fill(fill);
ee87a008 667 }
668 }
669 }
670
671}
672
673//________________________________________________________________________
674Int_t AliAnalysisTaskHJetEmbed::FindGeoMatchedJet(const AliEmcalJet* jet, const TClonesArray *jetArray, Double_t &dR)
675{
676 dR = 999;
677 if(!jetArray) return -1;
678
679 Int_t index = -1;
680 Int_t nJets = jetArray->GetEntries();
681 Double_t dRMax = 1;
682 for(Int_t ij=0; ij<nJets; ij++)
683 {
684 AliEmcalJet* jetTmp = dynamic_cast<AliEmcalJet*>(jetArray->At(ij));
aaea4d43 685 if(!jetTmp) continue;
ee87a008 686 if(TMath::Abs(jetTmp->Eta())>1) continue; // Generous eta cut
687 Double_t dPhi = GetDPhi(jet->Phi(),jetTmp->Phi());
688 Double_t dEta = jet->Eta()-jetTmp->Eta();
689 Double_t dRTmp = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
690 if(dRTmp<dRMax)
691 {
692 dRMax = dRTmp;
693 index = ij;
694 }
695 }
696 dR = dRMax;
697 return index;
698}
699
700//________________________________________________________________________
aaea4d43 701Int_t AliAnalysisTaskHJetEmbed::FindEnergyMatchedJet(const AliEmcalJet* jet, const TClonesArray *jetArray, Double_t &dR, Double_t &fraction)
ee87a008 702{
703 dR = 999;
aaea4d43 704 fraction=-1;
705 if(!jetArray || !jet) return -1;
ee87a008 706
707 Int_t index = -1;
708 Int_t nJets = jetArray->GetEntries();
aaea4d43 709 Double_t maxFrac = 0;
ee87a008 710 Int_t nJetC = (Int_t)jet->GetNumberOfConstituents();
aaea4d43 711 Double_t jetPt = jet->Pt();
ee87a008 712 for(Int_t ij=0; ij<nJets; ij++)
713 {
714 AliEmcalJet* jetTmp = dynamic_cast<AliEmcalJet*>(jetArray->At(ij));
aaea4d43 715 if(!jetTmp) continue;
ee87a008 716 if(TMath::Abs(jetTmp->Eta())>1) continue; // Generous eta cut
aaea4d43 717 if(GetJetDistance(jet,jetTmp)>1) continue;
718
ee87a008 719 Int_t nc = (Int_t)jetTmp->GetNumberOfConstituents();
720 Double_t sumPt = 0;
721 for(Int_t ic=0; ic<nc; ic++)
722 {
ee87a008 723 for(Int_t ijc=0; ijc<nJetC; ijc++)
724 {
aaea4d43 725 if(jetTmp->TrackAt(ic)==jet->TrackAt(ijc))
726 {
727 AliVParticle *part = (AliVParticle*)jet->TrackAt(ijc,fTrackArray);
728 sumPt += part->Pt();
729 }
ee87a008 730 }
731 }
732 Double_t frac = sumPt/jetPt;
aaea4d43 733 if(frac>maxFrac)
ee87a008 734 {
aaea4d43 735 maxFrac = frac;
ee87a008 736 index = ij;
737 }
738 }
aaea4d43 739 fraction = maxFrac;
ee87a008 740
741 if(index>0)
742 {
743 AliEmcalJet* jetTmp = dynamic_cast<AliEmcalJet*>(jetArray->At(index));
aaea4d43 744 if(jetTmp)
745 dR = GetJetDistance(jet,jetTmp);
ee87a008 746 }
747 return index;
748}
749
750//________________________________________________________________________
751Double_t AliAnalysisTaskHJetEmbed::CalculateDPhi(const Double_t phi1, const Double_t phi2)
752{
753 Double_t dPhi = phi1-phi2;
754 if(dPhi>2*pi) dPhi -= 2*pi;
755 if(dPhi<-2*pi) dPhi += 2*pi;
756 if(dPhi<-0.5*pi) dPhi += 2*pi;
757 if(dPhi>1.5*pi) dPhi -= 2*pi;
758 return dPhi;
759}
760
761//________________________________________________________________________
762Double_t AliAnalysisTaskHJetEmbed::GetDPhi(const Double_t phi1, const Double_t phi2)
763{
764 Double_t dPhi = TMath::Abs(phi1-phi2);
765 if(dPhi>2*pi) dPhi -= 2*pi;
766 if(dPhi>pi) dPhi = 2*pi - dPhi;
767 return dPhi;
768}
769
aaea4d43 770//________________________________________________________________________
771Double_t AliAnalysisTaskHJetEmbed::GetJetDistance(const AliEmcalJet *jet1, const AliEmcalJet* jet2)
772{
773 Double_t dPhi = GetDPhi(jet1->Phi(),jet2->Phi());
774 Double_t dEta = jet1->Eta()-jet2->Eta();
775 return TMath::Sqrt(dPhi*dPhi+dEta*dEta);
776}
777
ee87a008 778//________________________________________________________________________
779void AliAnalysisTaskHJetEmbed::RunQA()
780{
781 if(!fPLJetArray)
782 {
783 AliWarning(Form("Particle-level jet array is not available: %s\n",fPLJetArrName.Data()));
784 }
785 else
786 {
787 Int_t nPLJets = fPLJetArray->GetEntries();
788 for(Int_t ij=0; ij<nPLJets; ij++)
789 {
790 AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fPLJetArray->At(ij));
aaea4d43 791 if(!jet) continue;
ee87a008 792 if(!IsGoodJet(jet)) continue; // eta cut
793 Double_t jetPt = jet->Pt();
2942f542 794 Double_t fill[] = {jetPt, fCentrality, static_cast<Double_t>(fPtHardBin)};
ee87a008 795 fhPLJetPtVsCent[fTriggerType]->Fill(fill);
796 AliDebug(5, Form("PL jet %d has (pt,eta,phi) = (%2.2f,%2.2f,%2.2f)",ij,jetPt,jet->Eta(),jet->Phi()));
797 }
798 }
799
800 if(!fDLJetArray)
801 {
aaea4d43 802 AliWarning(Form("Detector-level jet array is not available: %s\n",fDLJetArrName.Data()));
ee87a008 803 }
804 else
805 {
806 Int_t nDLJets = fDLJetArray->GetEntries();
807 for(Int_t ij=0; ij<nDLJets; ij++)
808 {
809 AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fDLJetArray->At(ij));
aaea4d43 810 if(!jet) continue;
ee87a008 811 if(!IsGoodJet(jet)) continue; // eta cut
812 Double_t jetPt = jet->Pt();
2942f542 813 Double_t fill[] = {jetPt, fCentrality, static_cast<Double_t>(fPtHardBin)};
ee87a008 814 fhDLJetPtVsCent[fTriggerType]->Fill(fill);
815 AliDebug(5, Form("DL jet %d has pt = %2.2f",ij,jetPt));
816 }
817 }
818}
819
820
821//________________________________________________________________________
822Bool_t AliAnalysisTaskHJetEmbed::AcceptTrack(const AliVParticle *track)
823{
824 if(track->Pt()<fMinTrkPt || track->Pt()>fMaxTrkPt) return kFALSE;
825 if(track->Eta()<fMinTrkEta || track->Eta()>fMaxTrkEta) return kFALSE;
826 if(track->Phi()<fMinTrkPhi || track->Phi()>fMaxTrkPhi) return kFALSE;
827 return kTRUE;
828}
829
830//________________________________________________________________________
831Bool_t AliAnalysisTaskHJetEmbed::IsGoodJet(const AliEmcalJet* jet)
832{
833 Double_t etaCut = (0.9-fRadius>0.5)?0.5:0.9-fRadius;
834 if(TMath::Abs(jet->Eta())>etaCut) return kFALSE;
835 return kTRUE;
836}
837
838//
839//________________________________________________________________________
840//
841void AliAnalysisTaskHJetEmbed::PrintConfig()
842{
843 const char *decision[2] = {"no","yes"};
3ec095a2 844 const char *TTtype[2] = {"Single inclusive","Leading"};
ee87a008 845 printf("\n\n===== h-jet analysis configuration =====\n");
846 printf("Input event type: %s - %s\n",fCollisionSystem.Data(),fPeriod.Data());
847 printf("Track pt range: %2.2f < pt < %2.2f\n",fMinTrkPt, fMaxTrkPt);
848 printf("Track eta range: %2.1f < eta < %2.1f\n",fMinTrkEta, fMaxTrkEta);
849 printf("Track phi range: %2.0f < phi < %2.0f\n",fMinTrkPhi*TMath::RadToDeg(),fMaxTrkPhi*TMath::RadToDeg());
aaea4d43 850 printf("TT type: %s\n", TTtype[fTTtype]);
851 for(Int_t i=0; i<kNTT; i++)
852 printf("TT range %d: %2.0f < pt < %2.0f\n", i+1, fMinTTPt[i], fMaxTTPt[i]);
ee87a008 853 printf("Run QA: %s\n",decision[fRunQA]);
aaea4d43 854 printf("Run particle level: %s\n",decision[fRunPL]);
855 printf("Run detector level: %s\n",decision[fRunDL]);
ee87a008 856 printf("Run h+jet: %s\n",decision[fRunHJet]);
857 printf("Run matching: %s\n",decision[fRunMatch]);
858 printf("=======================================\n\n");
859}
860
861//________________________________________________________________________
862Double_t AliAnalysisTaskHJetEmbed::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)
863{
864 return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/(jetPx*jetPx+jetPy*jetPy+jetPz*jetPz);
865}
866
867//________________________________________________________________________
868void AliAnalysisTaskHJetEmbed::Terminate(Option_t *)
869{
870 // Called once at the end of the query
871}