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