Changes to compile with Root6 on macosx64
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskHJetEmbed.cxx
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 #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
45 ClassImp(AliAnalysisTaskHJetEmbed)
46
47 const Double_t pi = TMath::Pi();
48 const Double_t areaCut[4] = {0.1, 0.23, 0.4, 0.63};
49
50 //________________________________________________________________________
51 AliAnalysisTaskHJetEmbed::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), 
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 //________________________________________________________________________
107 AliAnalysisTaskHJetEmbed::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), 
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 //________________________________________________________________________
159 AliAnalysisTaskHJetEmbed::~AliAnalysisTaskHJetEmbed()
160 {
161   //Destructor
162   if(fRho)         delete fRho;
163   if(fOutputList) { fOutputList->Delete(); delete fOutputList;}
164 }
165
166 //________________________________________________________________________
167 void 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
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);
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
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 //________________________________________________________________________
305 void AliAnalysisTaskHJetEmbed::UserExec(Option_t *) 
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
482   fhCentrality[fTriggerType]->Fill(fCentrality);
483   fhRhoVsCent[fTriggerType]->Fill(fCentrality,fRhoValue);
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 //________________________________________________________________________
499 void AliAnalysisTaskHJetEmbed::RunHJet(const Double_t minPt, const Double_t maxPt)
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 //________________________________________________________________________
751 Double_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 //________________________________________________________________________
762 Double_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
770 //________________________________________________________________________
771 Double_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
778 //________________________________________________________________________
779 void 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));
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 //________________________________________________________________________
822 Bool_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 //________________________________________________________________________
831 Bool_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 //
841 void AliAnalysisTaskHJetEmbed::PrintConfig()
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 //________________________________________________________________________
868 void AliAnalysisTaskHJetEmbed::Terminate(Option_t *) 
869 {
870   // Called once at the end of the query
871 }