AliPhysics  251aa1e (251aa1e)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskHJetEmbed.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 <TParameter.h>
18 
19 #include "AliAODEvent.h"
20 #include "AliAODInputHandler.h"
21 #include "AliAnalysisManager.h"
22 #include "AliAnalysisTask.h"
23 #include "AliCentrality.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 #include "AliNamedString.h"
37 
38 #include "AliEmcalJet.h"
39 
40 #include <iostream>
41 using std::cout;
42 using std::cerr;
43 using std::endl;
44 
46 
47 const Double_t pi = TMath::Pi();
48 //const Double_t areaCut[4] = {0.1, 0.23, 0.4, 0.63};
49 
50 //________________________________________________________________________
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),
58  fTTtype(0),
59  fRadius(0.4), fJetArrName(""), fPLJetArrName(""), fDLJetArrName(""),
60  fJetArray(0x0), fPLJetArray(0x0), fDLJetArray(0x0),
61  fRhoName(""), fRho(0x0), fRhoValue(0),
62  fPtHardBinName(0x0), fPtHardBin(-1), fRandom(0x0),
63  fRunQA(kTRUE), fRunHJet(kTRUE), fRunMatch(kTRUE),
64  fRunPL(kFALSE), fRunDL(kTRUE),
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  }
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  }
104 }
105 
106 //________________________________________________________________________
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),
114  fTTtype(0),
115  fRadius(0.4), fJetArrName(""), fPLJetArrName(""), fDLJetArrName(""),
116  fJetArray(0x0), fPLJetArray(0x0), fDLJetArray(0x0),
117  fRhoName(""), fRho(0x0), fRhoValue(0),
118  fPtHardBinName(0x0), fPtHardBin(-1), fRandom(0x0),
119  fRunQA(kTRUE), fRunHJet(kTRUE), fRunMatch(kTRUE),
120  fRunPL(kFALSE), fRunDL(kTRUE),
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  }
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  }
157 }
158 //________________________________________________________________________
160 {
161  //Destructor
162  if(fRho) delete fRho;
163  if(fOutputList) { fOutputList->Delete(); delete fOutputList;}
164 }
165 
166 //________________________________________________________________________
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
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;
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};
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 
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};
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);
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 
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 
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);
279  fOutputList->Add(fhJetPhiGeoMatch[i]);
280 
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);
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 
298  fRandom = new TRandom3();
299 
300  PrintConfig();
301  PostData(1, fOutputList);
302 }
303 
304 //________________________________________________________________________
306 {
307  // Main loop, called for each event.
308 
309  AliDebug(5,"Entering UserExec");
310  fTriggerType = -1;
311  fEvent = InputEvent();
312  if (!fEvent)
313  {
314  AliError("Input event not available");
315  return;
316  }
317  AliDebug(5,"Got the input event");
318  if(!fPtHardBinName)
319  {
320  // Get embedded pt hard bin number
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();
332  fhEventStat->Fill(0.5);
333  UInt_t trigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
334 
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 
403  // Get MC particle array
404  if (fRunPL && !fMCParticleArray)
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
450  if (fRunPL && !fPLJetArray && !fPLJetArrName.IsNull())
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
466  if (fRunDL && !fDLJetArray && !fDLJetArrName.IsNull())
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 
484  fhPtHardBins->Fill(fPtHardBin);
485 
486  if(fRunQA) RunQA();
487  if(fRunHJet)
488  {
489  for(Int_t i=0; i<kNTT; i++)
490  RunHJet(fMinTTPt[i],fMaxTTPt[i]);
491  }
492 
493  PostData(1, fOutputList);
494  return;
495 }
496 
497 
498 //________________________________________________________________________
500 {
501  TArrayI arr;
502  Int_t counter = 0;
503  Int_t indexPL = -1;
504 
505  if(fRunPL)
506  {
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++)
513  {
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
527  {
528  if(maxPLPt<pt)
529  {
530  maxPLPt = pt;
531  indexPL = iPart;
532  }
533  }
534  }
535  arr.Set(counter);
536  if(fTTtype==0)
537  {
538  if(counter==0) indexPL = -1;
539  else if(counter==1) indexPL = arr.At(0);
540  else
541  {
542  Double_t pro = fRandom->Uniform() * counter;
543  indexPL = arr.At(TMath::FloorNint(pro));
544  }
545  }
546  arr.Reset();
547  }
548 
549 
550  // Find trigger track on detector level and after embedding
551  const Int_t Ntracks = fTrackArray->GetEntries();
552  Double_t maxDLPt = 0;
553  Int_t indexDL = -1;
554  arr.Set(Ntracks);
555  counter = 0;
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;
561  if(t->GetLabel()!=0)
562  {
563  //cout<<iTracks<<" "<<t->Pt()<<" "<<t->GetLabel()<<endl;
564  Double_t pt = t->Pt();
565  if(fTTtype==0)
566  {
567  if (pt<maxPt && pt>=minPt)
568  {
569  arr.AddAt(iTracks,counter);
570  counter++;
571  }
572  }
573  else if(fTTtype==1)
574  {
575  if(maxDLPt<pt)
576  {
577  maxDLPt = pt;
578  indexDL = iTracks;
579  }
580  }
581  }
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);
588  else
589  {
590  Double_t pro = fRandom->Uniform() * counter;
591  indexDL = arr.At(TMath::FloorNint(pro));
592  }
593  }
594  arr.Reset();
595 
596  AliDebug(2,Form("TT indices: PL=%d, DL=%d\n",indexPL,indexDL));
597 
598  // Run h+jet
599  if(fRunPL) FillHJetCor(fMCParticleArray, indexPL, fPLJetArray, fhPLTT[fTriggerType], fhPLHJet[fTriggerType], kFALSE);
600  if(fRunDL) FillHJetCor(fTrackArray, indexDL, fDLJetArray, fhDLTT[fTriggerType], fhDLHJet[fTriggerType], kFALSE);
601  FillHJetCor(fTrackArray, indexDL, fJetArray, fhTTPt[fTriggerType], fhHJet[fTriggerType], kTRUE);
602 
603  if(fRunMatch) RunMatch(fTrackArray, indexDL);
604 }
605 
606 //________________________________________________________________________
607 void AliAnalysisTaskHJetEmbed::FillHJetCor(const TClonesArray *tracks, const Int_t leadingIndex, const TClonesArray *jetArray, THnSparse *hTT, THnSparse *hn, Bool_t isBkg)
608 {
609  if(leadingIndex<0) return;
610 
611  AliVParticle *tt = (AliVParticle*) tracks->At(leadingIndex);
612  Double_t triggerPt = tt->Pt();
613  Double_t fill1[] = {triggerPt, fCentrality, static_cast<Double_t>(fPtHardBin) };
614  hTT->Fill(fill1);
615  AliDebug(2,Form("Found a trigger with pt = %2.2f",triggerPt));
616 
617  Double_t triggerPhi = tt->Phi();
618  if(triggerPhi<0) triggerPhi += 2*pi;
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));
623  if(!jet) continue;
624  if(!IsGoodJet(jet)) continue; // eta cut
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);
629  Double_t fill[] = {triggerPt,jetPt-jetArea*fRhoValue,dPhi,jetArea,fCentrality,static_cast<Double_t>(fPtHardBin)};
630  if(!isBkg) fill[1] = jetPt;
631  AliDebug(10,"Fill the histograms");
632  hn->Fill(fill);
633  }
634 }
635 
636 //________________________________________________________________________
637 void AliAnalysisTaskHJetEmbed::RunMatch(const TClonesArray *tracks, const Int_t leadingIndex)
638 {
639  if(leadingIndex<0) return;
640 
641  if(!fDLJetArray || !fJetArray)
642  {
643  AliWarning("Jet array is not available.");
644  return;
645  }
646 
647  AliVParticle *tt = (AliVParticle*) tracks->At(leadingIndex);
648  Double_t dR = 999, fraction = -1;
649  Int_t nJets = fDLJetArray->GetEntries();
650  for(Int_t ij=0; ij<nJets; ij++)
651  {
652  AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fDLJetArray->At(ij));
653  if(!jet) continue;
654  if(!IsGoodJet(jet)) continue; // eta cut
655  Double_t jetPt = jet->Pt();
656  if(jetPt<10) continue;
657 
658  // energy matching
659  Int_t mthJetIndexEn = FindEnergyMatchedJet(jet,fJetArray,dR,fraction);
660  if(mthJetIndexEn>-1 && fraction>0.5)
661  {
662  AliEmcalJet* jetMthEn = dynamic_cast<AliEmcalJet*>(fJetArray->At(mthJetIndexEn));
663  if(jetMthEn)
664  {
665  Double_t fill[] = {tt->Pt(),jetPt,CalculateDPhi(tt->Phi(),jet->Phi()),CalculateDPhi(jetMthEn->Phi(),jet->Phi()),dR,fCentrality,static_cast<Double_t>(fPtHardBin)};
666  fhJetPhiEnMatch[fTriggerType]->Fill(fill);
667  }
668  }
669  }
670 
671 }
672 
673 //________________________________________________________________________
674 Int_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));
685  if(!jetTmp) continue;
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 //________________________________________________________________________
701 Int_t AliAnalysisTaskHJetEmbed::FindEnergyMatchedJet(const AliEmcalJet* jet, const TClonesArray *jetArray, Double_t &dR, Double_t &fraction)
702 {
703  dR = 999;
704  fraction=-1;
705  if(!jetArray || !jet) return -1;
706 
707  Int_t index = -1;
708  Int_t nJets = jetArray->GetEntries();
709  Double_t maxFrac = 0;
710  Int_t nJetC = (Int_t)jet->GetNumberOfConstituents();
711  Double_t jetPt = jet->Pt();
712  for(Int_t ij=0; ij<nJets; ij++)
713  {
714  AliEmcalJet* jetTmp = dynamic_cast<AliEmcalJet*>(jetArray->At(ij));
715  if(!jetTmp) continue;
716  if(TMath::Abs(jetTmp->Eta())>1) continue; // Generous eta cut
717  if(GetJetDistance(jet,jetTmp)>1) continue;
718 
719  Int_t nc = (Int_t)jetTmp->GetNumberOfConstituents();
720  Double_t sumPt = 0;
721  for(Int_t ic=0; ic<nc; ic++)
722  {
723  for(Int_t ijc=0; ijc<nJetC; ijc++)
724  {
725  if(jetTmp->TrackAt(ic)==jet->TrackAt(ijc))
726  {
727  AliVParticle *part = (AliVParticle*)jet->TrackAt(ijc,fTrackArray);
728  sumPt += part->Pt();
729  }
730  }
731  }
732  Double_t frac = sumPt/jetPt;
733  if(frac>maxFrac)
734  {
735  maxFrac = frac;
736  index = ij;
737  }
738  }
739  fraction = maxFrac;
740 
741  if(index>0)
742  {
743  AliEmcalJet* jetTmp = dynamic_cast<AliEmcalJet*>(jetArray->At(index));
744  if(jetTmp)
745  dR = GetJetDistance(jet,jetTmp);
746  }
747  return index;
748 }
749 
750 //________________________________________________________________________
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 //________________________________________________________________________
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 
770 //________________________________________________________________________
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 
778 //________________________________________________________________________
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));
791  if(!jet) continue;
792  if(!IsGoodJet(jet)) continue; // eta cut
793  Double_t jetPt = jet->Pt();
794  Double_t fill[] = {jetPt, fCentrality, static_cast<Double_t>(fPtHardBin)};
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  {
802  AliWarning(Form("Detector-level jet array is not available: %s\n",fDLJetArrName.Data()));
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));
810  if(!jet) continue;
811  if(!IsGoodJet(jet)) continue; // eta cut
812  Double_t jetPt = jet->Pt();
813  Double_t fill[] = {jetPt, fCentrality, static_cast<Double_t>(fPtHardBin)};
814  fhDLJetPtVsCent[fTriggerType]->Fill(fill);
815  AliDebug(5, Form("DL jet %d has pt = %2.2f",ij,jetPt));
816  }
817  }
818 }
819 
820 
821 //________________________________________________________________________
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 //________________________________________________________________________
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 //
842 {
843  const char *decision[2] = {"no","yes"};
844  const char *TTtype[2] = {"Single inclusive","Leading"};
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());
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]);
853  printf("Run QA: %s\n",decision[fRunQA]);
854  printf("Run particle level: %s\n",decision[fRunPL]);
855  printf("Run detector level: %s\n",decision[fRunDL]);
856  printf("Run h+jet: %s\n",decision[fRunHJet]);
857  printf("Run matching: %s\n",decision[fRunMatch]);
858  printf("=======================================\n\n");
859 }
860 
861 //________________________________________________________________________
862 Double_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 //________________________________________________________________________
869 {
870  // Called once at the end of the query
871 }
Bool_t AcceptTrack(const AliVParticle *track)
Double_t CalculateDPhi(const Double_t phi1, const Double_t phi2)
Double_t Area() const
Definition: AliEmcalJet.h:117
Double_t fMinTrkPt
Index of the trigger track in the event.
double Double_t
Definition: External.C:58
TClonesArray * fDLJetArray
Array of the embedded PYTHIA jet array on particle level.
Definition: External.C:236
Double_t GetDPhi(const Double_t phi1, const Double_t phi2)
Double_t Eta() const
Definition: AliEmcalJet.h:108
TString fileName
Double_t Phi() const
Definition: AliEmcalJet.h:104
Double_t GetJetDistance(const AliEmcalJet *jet1, const AliEmcalJet *jet2)
centrality
TClonesArray * fPLJetArray
Array of the found jets.
TString fRhoName
Array of the embedded PYTHIA jet array on detector level.
Int_t fCollisionSystem
UShort_t GetNumberOfConstituents() const
Definition: AliEmcalJet.h:127
AliNamedString * fPtHardBinName
Value of the rho parameter.
int Int_t
Definition: External.C:63
Double_t fMaxVtxZ
V0M for current event.
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
void RunHJet(const Double_t minPt, const Double_t maxPt)
Int_t fPtHardBin
Pt hard bin param.
Int_t FindGeoMatchedJet(const AliEmcalJet *jet, const TClonesArray *jetArray, Double_t &dR)
Int_t FindEnergyMatchedJet(const AliEmcalJet *jet, const TClonesArray *jetArray, Double_t &dR, Double_t &fraction)
Double_t Pt() const
Definition: AliEmcalJet.h:96
Double_t fRhoValue
Rho parameter.
THnSparse * fhJetPhiGeoMatch[kNTrig]
Short_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:147
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_t fCentrality
Trigger type of the event.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
const char Option_t
Definition: External.C:48
const Double_t pi
bool Bool_t
Definition: External.C:53
void FillHJetCor(const TClonesArray *tracks, const Int_t leadingIndex, const TClonesArray *jetArray, THnSparse *hTT, THnSparse *hn, Bool_t isBkg=kFALSE)
void RunMatch(const TClonesArray *tracks, const Int_t leadingIndex)
Bool_t IsGoodJet(const AliEmcalJet *jet)
ClassImp(AliAnalysisTaskHJetEmbed) const Double_t pi
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65