AliPhysics  2c6b7ad (2c6b7ad)
AliAnalysisTaskJetCoreEmcal.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2016, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 #include <TClonesArray.h>
17 #include <TList.h>
18 #include "TChain.h"
19 #include "TTree.h"
20 #include "TMath.h"
21 #include "TH1F.h"
22 #include "TH2F.h"
23 #include "TH3F.h"
24 #include "THnSparse.h"
25 #include "TCanvas.h"
26 #include "TRandom3.h"
27 #include "AliLog.h"
28 
29 #include <AliAnalysisManager.h>
30 #include <AliVEventHandler.h>
31 #include <AliVEvent.h>
32 #include <AliVCluster.h>
33 #include <AliVParticle.h>
34 #include <AliLog.h>
35 #include <AliAnalysisDataSlot.h>
36 #include <AliAnalysisDataContainer.h>
37 
38 #include "AliTLorentzVector.h"
39 #include "AliEmcalJet.h"
40 #include "AliRhoParameter.h"
41 #include "AliJetContainer.h"
42 #include "AliParticleContainer.h"
43 #include "AliClusterContainer.h"
44 #include "AliEmcalPythiaInfo.h"
45 
47 
51 
57  fHistManager(),
58  fJetShapeType(AliAnalysisTaskJetCoreEmcal::kData),
59  fCentMin(0.),
60  fCentMax(100.),
61  fTTLowRef(8.),
62  fTTUpRef(9.),
63  fTTLowSig(20.),
64  fTTUpSig(50.),
65  fNRPBins(50),
66  fFrac(0.8),
67  fJetHadronDeltaPhi(0.6),
68  fMinFractionSharedPt(0.5),
69  fMinEmbJetPt(15.),
70  fJetContName(""),
71  fJetContTrueName(""),
72  fJetContPartName(""),
73  fFillTrackHistograms(kTRUE),
74  fFillJetHistograms(kTRUE),
75  fFillRecoilTHnSparse(kTRUE),
76  fFillInclusiveTree(kFALSE),
77  fFillRecoilTree(kFALSE),
78  fPtHardBin(0.),
79  fRejectionFactorInclusiveJets(1),
80  fRandom(0),
81  fHistEvtSelection(0x0),
82  fHJetSpec(0x0),
83  fh1TrigRef(0x0),
84  fh1TrigSig(0x0),
85  fh2Ntriggers(0x0),
86  fhDphiPtSig(0x0),
87  fhDphiPtRef(0x0),
88  fhPtDetPart(0x0),
89  fhPtHybrDet(0x0),
90  fhPtHybrPart(0x0),
91  fhPtHybrPartCor(0x0),
92  fhPhiHybrPartCor(0x0),
93  fhPtDet(0x0),
94  fhPtDetMatchedToPart(0x0),
95  fhResidual(0x0),
96  fhPtResidual(0x0),
97  fhPhiResidual(0x0),
98  fhPhiPhiResidual(0x0),
99  fhPtDetPartRecoil(0x0),
100  fhPtHybrDetRecoil(0x0),
101  fhPtHybrPartRecoil(0x0),
102  fhPtHybrPartCorRecoil(0x0),
103  fhPtDetRecoil(0x0),
104  fhPtDetMatchedToPartRecoil(0x0),
105  fhResidualRecoil(0x0),
106  fhPtResidualRecoil(0x0),
107  fhDphiResidualRecoil(0x0),
108  fhDphiphiResidualRecoil(0x0),
109  fhTTPtDetMatchedToPart(0x0),
110  fhTTPhiDetMatchedToPart(0x0),
111  fhDPhiHybrPartCorRecoil(0x0),
112  fTreeEmbInclusive(0x0),
113  fTreeEmbRecoil(0x0)
114 {
116 
117  DefineOutput(1, TList::Class());
118  //if(fJetShapeType == AliAnalysisTaskJetCoreEmcal::kDetEmbPart && fFillInclusiveTree)
119  DefineOutput(2, TTree::Class());
120 // if(fJetShapeType == AliAnalysisTaskJetCoreEmcal::kDetEmbPart && fFillRecoilTree)
121  DefineOutput(3, TTree::Class());
122 }
123 
130  AliAnalysisTaskEmcalJet(name, kTRUE),
131  fHistManager(name),
133  fCentMin(0.),
134  fCentMax(100.),
135  fTTLowRef(8),
136  fTTUpRef(9.),
137  fTTLowSig(20.),
138  fTTUpSig(50.),
139  fNRPBins(50),
140  fFrac(0.8),
141  fJetHadronDeltaPhi(0.6),
143  fMinEmbJetPt(15.),
144  fJetContName(""),
145  fJetContTrueName(""),
146  fJetContPartName(""),
147  fFillTrackHistograms(kTRUE),
148  fFillJetHistograms(kTRUE),
149  fFillRecoilTHnSparse(kTRUE),
150  fFillInclusiveTree(kFALSE),
151  fFillRecoilTree(kFALSE),
152  fPtHardBin(0.),
154  fRandom(0),
155  fHistEvtSelection(0x0),
156  fHJetSpec(0x0),
157  fh1TrigRef(0x0),
158  fh1TrigSig(0x0),
159  fh2Ntriggers(0x0),
160  fhDphiPtSig(0x0),
161  fhDphiPtRef(0x0),
162  fhPtDetPart(0x0),
163  fhPtHybrDet(0x0),
164  fhPtHybrPart(0x0),
165  fhPtHybrPartCor(0x0),
166  fhPhiHybrPartCor(0x0),
167  fhPtDet(0x0),
169  fhResidual(0x0),
170  fhPtResidual(0x0),
171  fhPhiResidual(0x0),
172  fhPhiPhiResidual(0x0),
173  fhPtDetPartRecoil(0x0),
174  fhPtHybrDetRecoil(0x0),
175  fhPtHybrPartRecoil(0x0),
177  fhPtDetRecoil(0x0),
179  fhResidualRecoil(0x0),
180  fhPtResidualRecoil(0x0),
186  fTreeEmbInclusive(0x0),
187  fTreeEmbRecoil(0x0)
188 {
190 
191  DefineOutput(1, TList::Class());
192  DefineOutput(2, TTree::Class());
193  DefineOutput(3, TTree::Class());
194 }
195 
200 {
201 }
202 
208 {
210 
213  //AllocateClusterHistograms();
214  //AllocateCellHistograms();
216 
217  TIter next(fHistManager.GetListOfHistograms());
218  TObject* obj = 0;
219  while ((obj = next())) {
220  fOutput->Add(obj);
221  }
222 
223  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
224  if(fJetShapeType == AliAnalysisTaskJetCoreEmcal::kDetEmbPart && fFillInclusiveTree) PostData(2, fTreeEmbInclusive); // Post data for ALL output slots > 0 here.
225  if(fJetShapeType == AliAnalysisTaskJetCoreEmcal::kDetEmbPart && fFillRecoilTree) PostData(3, fTreeEmbRecoil); // Post data for ALL output slots > 0 here.
226 }
227 
228 /*
229  * This function allocates the histograms for basic EMCal cluster QA.
230  * A set of histograms (energy, eta, phi, number of cluster) is allocated
231  * per each cluster container and per each centrality bin.
232  */
234 {
235  TString histname;
236  TString histtitle;
237  TString groupname;
238  AliClusterContainer* clusCont = 0;
239  TIter next(&fClusterCollArray);
240  while ((clusCont = static_cast<AliClusterContainer*>(next()))) {
241  groupname = clusCont->GetName();
242  // Protect against creating the histograms twice
243  if (fHistManager.FindObject(groupname)) {
244  AliWarning(TString::Format("%s: Found groupname %s in hist manager. The cluster containers will be filled into the same histograms.", GetName(), groupname.Data()));
245  continue;
246  }
247  fHistManager.CreateHistoGroup(groupname);
248  for (Int_t cent = 0; cent < fNcentBins; cent++) {
249  histname = TString::Format("%s/histClusterEnergy_%d", groupname.Data(), cent);
250  histtitle = TString::Format("%s;#it{E}_{cluster} (GeV);counts", histname.Data());
251  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
252 
253  histname = TString::Format("%s/histClusterEnergyExotic_%d", groupname.Data(), cent);
254  histtitle = TString::Format("%s;#it{E}_{cluster}^{exotic} (GeV);counts", histname.Data());
255  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
256 
257  histname = TString::Format("%s/histClusterNonLinCorrEnergy_%d", groupname.Data(), cent);
258  histtitle = TString::Format("%s;#it{E}_{cluster}^{non-lin.corr.} (GeV);counts", histname.Data());
259  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
260 
261  histname = TString::Format("%s/histClusterHadCorrEnergy_%d", groupname.Data(), cent);
262  histtitle = TString::Format("%s;#it{E}_{cluster}^{had.corr.} (GeV);counts", histname.Data());
263  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
264 
265  histname = TString::Format("%s/histClusterPhi_%d", groupname.Data(), cent);
266  histtitle = TString::Format("%s;#it{#phi}_{custer};counts", histname.Data());
267  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
268 
269  histname = TString::Format("%s/histClusterEta_%d", groupname.Data(), cent);
270  histtitle = TString::Format("%s;#it{#eta}_{custer};counts", histname.Data());
271  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
272 
273  histname = TString::Format("%s/histNClusters_%d", groupname.Data(), cent);
274  histtitle = TString::Format("%s;number of clusters;events", histname.Data());
275  if (fForceBeamType != kpp) {
276  fHistManager.CreateTH1(histname, histtitle, 500, 0, 3000);
277  }
278  else {
279  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
280  }
281  }
282  }
283 
284  histname = "fHistSumNClusters";
285  histtitle = TString::Format("%s;Sum of n clusters;events", histname.Data());
286  if (fForceBeamType != kpp) {
287  fHistManager.CreateTH1(histname, histtitle, 500, 0, 3000);
288  }
289  else {
290  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
291  }
292 }
293 
294 /*
295  * This function allocates the histograms for basic EMCal QA.
296  * One 2D histogram with the cell energy spectra and the number of cells
297  * per event is allocated per each centrality bin.
298  */
300 {
301  TString histname;
302  TString histtitle;
303  TString groupname(fCaloCellsName);
304 
305  fHistManager.CreateHistoGroup(groupname);
306  for (Int_t cent = 0; cent < fNcentBins; cent++) {
307  histname = TString::Format("%s/histCellEnergy_%d", groupname.Data(), cent);
308  histtitle = TString::Format("%s;#it{E}_{cell} (GeV);counts", histname.Data());
309  fHistManager.CreateTH1(histname, histtitle, 300, 0, 150);
310 
311  histname = TString::Format("%s/histNCells_%d", groupname.Data(), cent);
312  histtitle = TString::Format("%s;number of cells;events", histname.Data());
313  if (fForceBeamType != kpp) {
314  fHistManager.CreateTH1(histname, histtitle, 500, 0, 6000);
315  }
316  else {
317  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
318  }
319  }
320 }
321 
322 /*
323  * This function allocates the histograms for basic tracking QA.
324  * A set of histograms (pT, eta, phi, difference between kinematic properties
325  * at the vertex and at the EMCal surface, number of tracks) is allocated
326  * per each particle container and per each centrality bin.
327  */
329 {
330  TString histname;
331  TString histtitle;
332  TString groupname;
333  AliParticleContainer* partCont = 0;
334  TIter next(&fParticleCollArray);
335  while ((partCont = static_cast<AliParticleContainer*>(next()))) {
336  groupname = partCont->GetName();
337 
338  // Protect against creating the histograms twice
339  if (fHistManager.FindObject(groupname)) {
340  AliWarning(TString::Format("%s: Found groupname %s in hist manager. The track containers will be filled into the same histograms.", GetName(), groupname.Data()));
341  continue;
342  }
343  fHistManager.CreateHistoGroup(groupname);
344  for (Int_t cent = 0; cent < fNcentBins; cent++) {
345  histname = TString::Format("%s/histTrackPt_%d", groupname.Data(), cent);
346  histtitle = TString::Format("%s;#it{p}_{T,track} (GeV/#it{c});counts", histname.Data());
347  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
348 
349  histname = TString::Format("%s/histTrackPhi_%d", groupname.Data(), cent);
350  histtitle = TString::Format("%s;#it{#phi}_{track};counts", histname.Data());
351  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
352 
353  histname = TString::Format("%s/histTrackEta_%d", groupname.Data(), cent);
354  histtitle = TString::Format("%s;#it{#eta}_{track};counts", histname.Data());
355  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
356 
357  if (TClass(partCont->GetClassName()).InheritsFrom("AliVTrack")) {
358  histname = TString::Format("%s/fHistDeltaEtaPt_%d", groupname.Data(), cent);
359  histtitle = TString::Format("%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{#eta}_{track}^{vertex} - #it{#eta}_{track}^{EMCal};counts", histname.Data());
360  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, 50, -0.5, 0.5);
361 
362  histname = TString::Format("%s/fHistDeltaPhiPt_%d", groupname.Data(), cent);
363  histtitle = TString::Format("%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{#phi}_{track}^{vertex} - #it{#phi}_{track}^{EMCal};counts", histname.Data());
364  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, 200, -2, 2);
365 
366  histname = TString::Format("%s/fHistDeltaPtvsPt_%d", groupname.Data(), cent);
367  histtitle = TString::Format("%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{p}_{T,track}^{vertex} - #it{p}_{T,track}^{EMCal} (GeV/#it{c});counts", histname.Data());
368  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, fNbins / 2, -fMaxBinPt/2, fMaxBinPt/2);
369 
370  histname = TString::Format("%s/fHistEoverPvsP_%d", groupname.Data(), cent);
371  histtitle = TString::Format("%s;#it{P}_{track} (GeV/#it{c});#it{E}_{cluster} / #it{P}_{track} #it{c};counts", histname.Data());
372  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, fNbins / 2, 0, 4);
373  }
374 
375  histname = TString::Format("%s/histNTracks_%d", groupname.Data(), cent);
376  histtitle = TString::Format("%s;number of tracks;events", histname.Data());
377  if (fForceBeamType != kpp) {
378  fHistManager.CreateTH1(histname, histtitle, 500, 0, 5000);
379  }
380  else {
381  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
382  }
383  }
384  }
385 
386  histname = "fHistSumNTracks";
387  histtitle = TString::Format("%s;Sum of n tracks;events", histname.Data());
388  if (fForceBeamType != kpp) {
389  fHistManager.CreateTH1(histname, histtitle, 500, 0, 5000);
390  }
391  else {
392  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
393  }
394 }
395 
396 /*
397  * This function allocates the histograms for basic jet QA.
398  * A set of histograms (pT, eta, phi, area, number of jets, corrected pT) is allocated
399  * per each jet container and per each centrality bin.
400  */
402 {
403  TString histname;
404  TString histtitle;
405  TString groupname;
406  AliJetContainer* jetCont = 0;
407  TIter next(&fJetCollArray);
408  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
409  groupname = jetCont->GetName();
410  // Protect against creating the histograms twice
411  if (fHistManager.FindObject(groupname)) {
412  AliWarning(TString::Format("%s: Found groupname %s in hist manager. The jet containers will be filled into the same histograms.", GetName(), groupname.Data()));
413  continue;
414  }
415  fHistManager.CreateHistoGroup(groupname);
416  for (Int_t cent = 0; cent < fNcentBins; cent++) {
417  histname = TString::Format("%s/histJetPt_%d", groupname.Data(), cent);
418  histtitle = TString::Format("%s;#it{p}_{T,jet} (GeV/#it{c});counts", histname.Data());
419  fHistManager.CreateTH1(histname, histtitle, fNbins, fMinBinPt, fMaxBinPt);
420 
421  histname = TString::Format("%s/histJetArea_%d", groupname.Data(), cent);
422  histtitle = TString::Format("%s;#it{A}_{jet};counts", histname.Data());
423  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, 3);
424 
425  histname = TString::Format("%s/histJetPhi_%d", groupname.Data(), cent);
426  histtitle = TString::Format("%s;#it{#phi}_{jet};counts", histname.Data());
427  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
428 
429  histname = TString::Format("%s/histJetEta_%d", groupname.Data(), cent);
430  histtitle = TString::Format("%s;#it{#eta}_{jet};counts", histname.Data());
431  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
432 
433  histname = TString::Format("%s/histNJets_%d", groupname.Data(), cent);
434  histtitle = TString::Format("%s;number of jets;events", histname.Data());
435  if (fForceBeamType != kpp) {
436  fHistManager.CreateTH1(histname, histtitle, 500, 0, 500);
437  }
438  else {
439  fHistManager.CreateTH1(histname, histtitle, 100, 0, 100);
440  }
441 
442  if (!jetCont->GetRhoName().IsNull()) {
443  histname = TString::Format("%s/histJetCorrPt_%d", groupname.Data(), cent);
444  histtitle = TString::Format("%s;#it{p}_{T,jet}^{corr} (GeV/#it{c});counts", histname.Data());
445  fHistManager.CreateTH1(histname, histtitle, fNbins, -fMaxBinPt / 2, fMaxBinPt / 2);
446  }
447  }
448  }
449 }
450 
452 {
453 
454  Bool_t oldStatus = TH1::AddDirectoryStatus();
455  TH1::AddDirectory(kFALSE);
456 
457  // set seed
458  fRandom = new TRandom3(0);
459 
460  fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
461  fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
462  fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
463  fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
464  fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
465  fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
466  fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
467 
468  fh1TrigRef=new TH1D("Trig Ref","",10,0.,10);
469  fh1TrigSig=new TH1D("Trig Sig","",10,0.,10);
470  fh2Ntriggers=new TH2F("# of triggers","",100,0.,100.,50,0.,50.);
471 
473 
474  fOutput->Add(fh1TrigRef);
475  fOutput->Add(fh1TrigSig);
476  fOutput->Add(fh2Ntriggers);
477 
479  const Int_t dimSpec = 6;
480  const Int_t nBinsSpec[dimSpec] = {100,100, 280, 50,200, fNRPBins};
481  const Double_t lowBinSpec[dimSpec] = {0,0,-80, 0,-0.5*TMath::Pi(), 0};
482  const Double_t hiBinSpec[dimSpec] = {100,1, 200, 50,1.5*TMath::Pi(), static_cast<Double_t>(fNRPBins)};
483  fHJetSpec = new THnSparseF("fHJetSpec","Recoil jet spectrum",dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
484  }
485 
486  // comment out since I want finer binning in jet area, to make it easier
487  // to change selection on jet area (Leticia used 0.8*R^2*Pi whereas 0.6 is used
488  // for inclusive jets)
489  //change binning in jet area
490 // Double_t *xPt6 = new Double_t[7];
491 // xPt6[0] = 0.;
492 // xPt6[1]=0.07;
493 // xPt6[2]=0.2;
494 // xPt6[3]=0.4;
495 // xPt6[4]=0.6;
496 // xPt6[5]=0.8;
497 // xPt6[6]=1;
498 // fHJetSpec->SetBinEdges(1,xPt6);
499 // delete [] xPt6;
500 
501  fOutput->Add(fHJetSpec);
502 
503  // azimuthal correlation
504 
505  fhDphiPtSig = new TH2F("hDphiPtS","recoil #Delta #phi vs jet pT signal",100,-2,5,250,-50,200);
506  fhDphiPtSig->GetXaxis()->SetTitle("#Delta #phi");
507  fhDphiPtSig->GetYaxis()->SetTitle("p^{reco,ch}_{T,jet} (GeV/c)");
508  fhDphiPtRef = new TH2F("hDphiPtR","recoil #Delta #phi vs jet pT reference",100,-2,5,250,-50,200);
509  fhDphiPtRef->GetXaxis()->SetTitle("#Delta #phi");
510  fhDphiPtRef->GetYaxis()->SetTitle("p^{reco,ch}_{T,jet} (GeV/c)");
511 
512  fOutput->Add(fhDphiPtRef);
513  fOutput->Add(fhDphiPtSig);
514 
515 
516  fhPtHybrDet= new TH2F("hPtHybrDet","pT response Pb-Pb+PYTHIA vs PYTHIA",200,0,200,200,0,200);
517  fhPtHybrDet->GetXaxis()->SetTitle("p^{Pb-Pb+PYTHIA,ch}_{T} (GeV/c)");
518  fhPtHybrDet->GetYaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
519  fhPtDetPart= new TH2F("hPtDetPart","pT response PYTHIA vs part",200,0,200,200,0,200);
520  fhPtDetPart->GetXaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
521  fhPtDetPart->GetYaxis()->SetTitle("p^{true}_{T} (GeV/c)");
522  fhPtHybrPart = new TH2F("hPtHybrPart",Form("pT response Pb-Pb+PYTHIA vs part, min shared pT > %.0f",fMinFractionSharedPt*100),200,0,200,200,0,200);
523  fhPtHybrPart->GetXaxis()->SetTitle("p^{Pb-Pb+PYTHIA,ch}_{T} (GeV/c)");
524  fhPtHybrPart->GetYaxis()->SetTitle("p^{true}_{T} (GeV/c)");
525  fhPtHybrPartCor = new TH2F("hPtHybrPartCor",Form("pT response Pb-Pb+PYTHIA corrected vs part, min shared pT > %.0f",fMinFractionSharedPt*100),200,0,200,200,0,200);
526  fhPtHybrPartCor->GetXaxis()->SetTitle("p^{Pb-Pb+PYTHIA,ch}_{T} (GeV/c)");
527  fhPtHybrPartCor->GetYaxis()->SetTitle("p^{true}_{T} (GeV/c)");
528 
529  fhPhiHybrPartCor = new TH2F("hPhiHybrPartCor",Form("phi response Pb-Pb+PYTHIA corrected vs part, min shared pT > %.0f",fMinFractionSharedPt*100),200,-0.5*TMath::Pi(),1.5*TMath::Pi(),200,-0.5*TMath::Pi(),1.5*TMath::Pi());
530  fhPhiHybrPartCor->GetXaxis()->SetTitle("#phi^{Pb-Pb+PYTHIA,ch}");
531  fhPhiHybrPartCor->GetYaxis()->SetTitle("#phi^{true}");
532  fhResidual = new TH1F("hResidual","residual",50,-1,1);
533  fhResidual->GetXaxis()->SetTitle("p^{reco}_{T} - p^{true}_{T} / p^{true}_{T}");
534  fhPtResidual= new TH2F("hPtResidual","pT vs residual",200,0,200,50,-1,1);
535  fhPtResidual->GetXaxis()->SetTitle("p^{true}_{T} (GeV/c)");
536  fhPtResidual->GetYaxis()->SetTitle("p^{reco}_{T} - p^{true}_{T} / p^{true}_{T}");
537  fhPhiResidual = new TH1F("hPhiResidual","residual phi",50,-1,1);
538  fhPhiResidual->GetXaxis()->SetTitle("#phi^{reco} - #phi^{true} / #phi^{true}");
539  fhPhiPhiResidual= new TH2F("hPhiPhiResidual","pT vs residual",600,-3,3,200,-0.5*TMath::Pi(),1.5*TMath::Pi());
540  fhPhiPhiResidual->GetXaxis()->SetTitle("#phi^{true}");
541  fhPhiPhiResidual->GetYaxis()->SetTitle("#phi^{reco} - #phi^{true} / #phi^{true}");
542 
543  fhPtDet= new TH1F("hPtDet","pT detector level",200,0,200);
544  fhPtDet->GetXaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
545  fhPtDetMatchedToPart = new TH1F("hPtDetMatchedToPart","pT detector level matched to particle level jet",200,0,200);
546  fhPtDetMatchedToPart->GetXaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
547 
548  fhPtHybrDetRecoil= new TH2F("hPtHybrDetRecoil","pT response Pb-Pb+PYTHIA vs PYTHIA",200,0,200,200,0,200);
549  fhPtHybrDetRecoil->GetXaxis()->SetTitle("p^{Pb-Pb+PYTHIA,ch}_{T} (GeV/c)");
550  fhPtHybrDetRecoil->GetYaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
551  fhPtDetPartRecoil= new TH2F("hPtDetPartRecoil","pT response PYTHIA vs part",200,0,200,200,0,200);
552  fhPtDetPartRecoil->GetXaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
553  fhPtDetPartRecoil->GetYaxis()->SetTitle("p^{true}_{T} (GeV/c)");
554  fhPtHybrPartRecoil = new TH2F("hPtHybrPartRecoil",Form("pT response Pb-Pb+PYTHIA vs part, min shared pT > %.0f",fMinFractionSharedPt*100),200,0,200,200,0,200);
555  fhPtHybrPartRecoil->GetXaxis()->SetTitle("p^{Pb-Pb+PYTHIA,ch}_{T} (GeV/c)");
556  fhPtHybrPartRecoil->GetYaxis()->SetTitle("p^{true}_{T} (GeV/c)");
557  fhPtHybrPartCorRecoil = new TH2F("hPtHybrPartCorRecoil",Form("pT response Pb-Pb+PYTHIA corrected vs part, min shared pT > %.0f",fMinFractionSharedPt*100),200,0,200,200,0,200);
558  fhPtHybrPartCorRecoil->GetXaxis()->SetTitle("p^{Pb-Pb+PYTHIA,ch}_{T} (GeV/c)");
559  fhPtHybrPartCorRecoil->GetYaxis()->SetTitle("p^{true}_{T} (GeV/c)");
560 
561  fhDPhiHybrPartCorRecoil = new TH2F("hDPhiHybrPartCorRecoil",Form("#Delta#phi response Pb-Pb+PYTHIA corrected vs part, min shared pT > %.0f",fMinFractionSharedPt*100),200,-0.5*TMath::Pi(),1.5*TMath::Pi(),200,-0.5*TMath::Pi(),1.5*TMath::Pi());
562  fhDPhiHybrPartCorRecoil->GetXaxis()->SetTitle("#Delta#phi^{Pb-Pb+PYTHIA,ch}");
563  fhDPhiHybrPartCorRecoil->GetYaxis()->SetTitle("#Delta#phi^{true}");
564  fhResidualRecoil = new TH1F("hResidualRecoil","residual",50,-1,1);
565  fhResidualRecoil->GetXaxis()->SetTitle("p^{reco}_{T} - p^{true}_{T} / p^{true}_{T}");
566  fhPtResidualRecoil= new TH2F("hPtResidualRecoil","pT vs residual",200,0,200,50,-1,1);
567  fhPtResidualRecoil->GetXaxis()->SetTitle("p^{true}_{T} (GeV/c)");
568  fhPtResidualRecoil->GetYaxis()->SetTitle("p^{reco}_{T} - p^{true}_{T} / p^{true}_{T}");
569 
570  fhDphiResidualRecoil = new TH1F("hDphiResidualRecoil","residual",50,-1,1);
571  fhDphiResidualRecoil->GetXaxis()->SetTitle("#Delta#phi^{reco} - #Delta#phi^{true} / #Delta#phi^{true}");
572  fhDphiphiResidualRecoil= new TH2F("hDphiphiResidualRecoil","#Delta#phi vs residual in phi",400,-3,3,200,-0.5*TMath::Pi(),1.5*TMath::Pi());
573  fhDphiphiResidualRecoil->GetXaxis()->SetTitle("#Delta#phi^{true}");
574  fhDphiphiResidualRecoil->GetYaxis()->SetTitle("#Delta#phi^{reco} - #Delta#phi^{true} / #Delta#phi^{true}");
575  fhPtDetRecoil= new TH1F("hPtDetRecoil","pT detector level",200,0,200);
576  fhPtDetRecoil->GetXaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
577  fhPtDetMatchedToPartRecoil = new TH1F("hPtDetMatchedToPartRecoil","pT detector level matched to particle level jet",200,0,200);
578  fhPtDetMatchedToPartRecoil->GetXaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
579 
580  fhTTPtDetMatchedToPart = new TH2F("hTTPtDetMatchedToPart","trigger track pT response reco vs partice",140,0,70,140,0,70);
581  fhTTPtDetMatchedToPart->GetXaxis()->SetTitle("p^{TT,reco}_{T} (GeV/c)");
582  fhTTPtDetMatchedToPart->GetYaxis()->SetTitle("p^{TT,part}_{T} (GeV/c)");
583  fhTTPhiDetMatchedToPart = new TH2F("hTTPhiDetMatchedToPart","trigger track #phi response reco vs partice",200,-0.5*TMath::Pi(),1.5*TMath::Pi(),200,-0.5*TMath::Pi(),1.5*TMath::Pi());
584  fhTTPhiDetMatchedToPart->GetXaxis()->SetTitle("p^{TT,reco}_{T} (GeV/c)");
585  fhTTPhiDetMatchedToPart->GetYaxis()->SetTitle("p^{TT,part}_{T} (GeV/c)");
586 
587  fOutput->Add(fhPtHybrDet);
588  fOutput->Add(fhPtDetPart);
589  fOutput->Add(fhPtHybrPart);
590  fOutput->Add(fhPtHybrPartCor);
592  fOutput->Add(fhPtDet);
594  fOutput->Add(fhResidual);
595  fOutput->Add(fhPtResidual);
596  fOutput->Add(fhPhiResidual);
598 
604  fOutput->Add(fhPtDetRecoil);
610 
613 
614  TString varNamesInclusive[8]={"centrality","ptRawRec","areaRec","ptCorrRec","phiRec","ptPart","phiPart","binPtHard"};
615  TString varNamesRecoil[11]={"centrality","ptTT","ptRawRec","areaRec","ptCorrRec","phiRec","DPhiRec","ptPart","phiPart","DPhiPart","binPtHard"};
617  const char* nameEmbInclusive = GetOutputSlot(2)->GetContainer()->GetName();
618  fTreeEmbInclusive = new TTree(nameEmbInclusive, nameEmbInclusive);
619  for(Int_t ivar=0; ivar < 8; ivar++){
620  fTreeEmbInclusive->Branch(varNamesInclusive[ivar].Data(), &fTreeVarsInclusive[ivar], Form("%s/F", varNamesInclusive[ivar].Data()));
621  }
622  }
623 
625  const char* nameEmbRecoil= GetOutputSlot(3)->GetContainer()->GetName();
626  fTreeEmbRecoil = new TTree(nameEmbRecoil, nameEmbRecoil);
627  for(Int_t ivar=0; ivar < 11; ivar++){
628  fTreeEmbRecoil->Branch(varNamesRecoil[ivar].Data(), &fTreeVarsRecoil[ivar], Form("%s/F", varNamesRecoil[ivar].Data()));
629  }
630  }
631 
632  // =========== Switch on Sumw2 for all histos ===========
633  for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
634  TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
635  if (h1){
636  h1->Sumw2();
637  continue;
638  }
639  THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
640  if (hn){
641  hn->Sumw2();
642  }
643  }
644 
645  // add QA plots from fEventCuts
646  fEventCuts.AddQAplotsToList(fOutput);
647 
648  TH1::AddDirectory(oldStatus);
649 
650  PostData(1, fOutput);
653 }
654 
662 {
663 
664  fHistEvtSelection->Fill(1); // number of events before event selection
665  AliVEvent *ev = InputEvent();
666  if (!fEventCuts.AcceptEvent(ev)) {
667  fHistEvtSelection->Fill(2);
668  return kTRUE;
669  }
670 
673  //DoClusterLoop();
674  //DoCellLoop();
675  DoJetCoreLoop();
677 
678  return kTRUE;
679 }
680 
682 {
683  // Do jet core analysis and fill histograms.
684 
686  AliJetContainer *jetContPart = 0x0;
687  AliJetContainer *jetContTrue = 0x0;
689  jetContTrue = GetJetContainer(fJetContTrueName);
690  jetContPart = GetJetContainer(fJetContPartName);
691  }
692 
693  if(!jetCont ||
695  AliError(Form("jet container not found - check name %s",fJetContName.Data()));
696  TIter next(&fJetCollArray);
697  while ((jetCont = static_cast<AliJetContainer*>(next())))
698  AliError(Form("%s",jetCont->GetName()));
699  AliFatal("Exit...");
700  return;
701  }
702 
703  // centrality selection
704  if(fDebug) Printf("centrality: %f\n", fCent);
706  ((fCent>fCentMax) || (fCent<fCentMin))) {
707  fHistEvtSelection->Fill(4);
708  return;
709  }
710  fHistEvtSelection->Fill(0);
711 
712  // Background
713  // If rho exists get it, otherwise it is set to 0 and ptJet_corr = ptJet_raw
714  Double_t rho = 0;
715  if (jetCont->GetRhoParameter()) rho = jetCont->GetRhoVal();
716  if(fDebug) Printf("rho = %f",rho);
717 
718  // get MC particle container in case running embedding, to match
719  // reconstructed and MC-level trigger tracks
720  AliParticleContainer *partCont = 0x0;
722 // AliParticleContainer *partCont = 0x0;
723 // TIter next(&fParticleCollArray);
724 // while ((partCont = static_cast<AliParticleContainer*>(next()))) {
725 // TString groupname = partCont->GetName();
726 // Printf("particle name = %s",groupname.Data());
727 // }
728 
729  // Choose trigger track
730  Int_t nT=0;
731  TList ParticleList;
732  Double_t minT=0;
733  Double_t maxT=0;
734  Int_t number=0;
735  Double_t dice=fRandom->Uniform(0,1);
736  Bool_t isSignal = kFALSE;
737  if(dice>fFrac){
738  minT=fTTLowRef;
739  maxT=fTTUpRef;
740  }
741  if(dice<=fFrac){
742  isSignal = kTRUE;
743  minT=fTTLowSig;
744  maxT=fTTUpSig;
745  }
746  nT=SelectTrigger(&ParticleList,minT,maxT,number);
747  if(fDebug) Printf("%s class ---> n triggers between %f and %f = %i, index of trigger chosen = %i",dice>fFrac?"ref.":"sig.",minT,maxT,number,nT);
748  if(nT<0) return;
749 
750  if(dice>fFrac) fh1TrigRef->Fill(number);
751  if(dice<=fFrac)fh1TrigSig->Fill(number);
752 
753 
754  // particle loop -
755  for(Int_t tt=0;tt<ParticleList.GetEntries();tt++){
756  // histogram
757  //if(fHardest==0||fHardest==1){if(tt!=nT) continue;}
758  if(tt!=nT) continue;
759  AliVParticle *partback = (AliVParticle*)ParticleList.At(tt);
760  if(!partback) continue;
761  if(fDebug) Printf("trigger particle pt = %f \teta = %f \t phi = %f",partback->Pt(),partback->Eta(),partback->Phi());
762  // if(partback->Pt()<8) continue;
763 
764  fh2Ntriggers->Fill(fCent,partback->Pt());
765  Double_t phiBinT = RelativePhi(partback->Phi(),fEPV0);
766 
767  Double_t etabig=0;
768  Double_t ptbig=0;
769  Double_t areabig=0;
770  Double_t phibig=0.;
771  // Double_t areasmall=0;
772 
773  TString histname;
774  TString groupname;
775  groupname = jetCont->GetName();
776  UInt_t count = 0;
777  for(auto jetbig : jetCont->accepted()) {
778  if (!jetbig) continue;
779  count++;
780  ptbig = jetbig->Pt();
781  etabig = jetbig->Eta();
782  phibig = jetbig->Phi();
783  if(ptbig==0) continue;
784  Double_t phiBin = RelativePhi(phibig,fEPV0); //relative phi between jet and ev. plane
785  areabig = jetbig->Area();
786  Double_t ptcorr=ptbig-rho*areabig;
787  Double_t dphi=RelativePhi(partback->Phi(),phibig);
788  if(fDebug) Printf("jet properties...\n\teta = %f \t phi = %f \t pt = %f \t relativephi = %f\t area = %f\t rho = %f",etabig,phibig,ptbig,dphi,areabig,rho);
789 
790  // do azimuthal correlation analysis
791  // dPhi between -0.5 < dPhi < 1.5
792  Double_t dPhiShift=phibig-partback->Phi();
793  if(dPhiShift>2*TMath::Pi()) dPhiShift -= 2*TMath::Pi();
794  if(dPhiShift<-2*TMath::Pi()) dPhiShift += 2*TMath::Pi();
795  if(dPhiShift<-0.5*TMath::Pi()) dPhiShift += 2*TMath::Pi();
796  if(dPhiShift>1.5*TMath::Pi()) dPhiShift -= 2*TMath::Pi();
797  if(isSignal) fhDphiPtSig->Fill(dPhiShift,ptcorr);
798  else fhDphiPtRef->Fill(dPhiShift,ptcorr);
799 
800  // selection on relative phi
801  if(fJetHadronDeltaPhi>0. &&
802  TMath::Abs(dphi)<TMath::Pi()-fJetHadronDeltaPhi) continue;
803 
805  Float_t phitt=partback->Phi();
806  if(phitt<0)phitt+=TMath::Pi()*2.;
807  Int_t phiBintt = GetPhiBin(phitt-fEPV0);
808 
809  Double_t fillspec[] = {fCent,areabig,ptcorr,partback->Pt(),dPhiShift, static_cast<Double_t>(phiBintt)};
810  fHJetSpec->Fill(fillspec);
811  }
812 
814  //
815  // embedding for recoil jets
816  // get MC info
817  //
818  if(ptcorr<fMinEmbJetPt) continue;
819  Double_t ptTTMC = 0;
820  Double_t phiTTMC = 0;
821  Int_t TTmatched = 0;
822  for(auto partMC : partCont->accepted()) {
823  Int_t labtr = partback->GetLabel();
824  Int_t labpa = partMC->GetLabel();
825  if(labtr==labpa) {
826  ptTTMC = partMC->Pt();
827  phiTTMC = partMC->Phi();
828  TTmatched++;
829  break;
830  }
831  }
832  if(TTmatched!=1) continue;
833  Double_t ptTTreco = partback->Pt();
834  Double_t phiTTreco = partback->Phi();
835  if(fDebug) Printf("found corresponding truth-level particle, pt reco = %f pt MC = %f, phi reco = %f phi MC = %f",ptTTreco,ptTTMC,phiTTreco,phiTTMC);
836  fhTTPtDetMatchedToPart->Fill(ptTTreco,ptTTMC);
837  fhTTPhiDetMatchedToPart->Fill(phiTTreco,phiTTMC);
838 
839  auto jet2 = jetbig->ClosestJet();
840  if(!jet2) {
841  //Printf("jet 2 cant be found");
842  continue;}
843  Double_t ptJet2 = jet2->Pt();
844  fhPtDetRecoil->Fill(ptJet2);
845  auto jet3 = jet2->ClosestJet();
846  if(!jet3) {
847  //Printf("jet3 can't be found");
848  continue;
849  }
850  fhPtDetMatchedToPartRecoil->Fill(ptJet2);
851  Double_t ptJet3 = jet3->Pt();
852  Double_t phiJet3 = jet3->Phi();
853 
854 
855  Double_t dPhiPart=phiJet3-phiTTMC;
856  if(dPhiPart>2*TMath::Pi()) dPhiPart -= 2*TMath::Pi();
857  if(dPhiPart<-2*TMath::Pi()) dPhiPart += 2*TMath::Pi();
858  if(dPhiPart<-0.5*TMath::Pi()) dPhiPart += 2*TMath::Pi();
859  if(dPhiPart>1.5*TMath::Pi()) dPhiPart -= 2*TMath::Pi();
860 
861  if(fDebug) Printf("--- recoil - jet pt = jet hybrid pt = %f\t jet matched det pt = %f\t jet matched particle level pt = %f\t\n\tjet reco phi = %f\t jet particle phi = %f",ptbig,ptJet2,ptJet3,phibig,phiJet3);
862 
863  fhPtDetPartRecoil->Fill(ptJet2,ptJet3);
864  Double_t fraction = jetCont->GetFractionSharedPt(jetbig);
865  if(fraction < fMinFractionSharedPt) continue;
866 
867  Double_t residual = (ptcorr - ptJet3) / ptJet3;
868  Double_t residualDphi = (dPhiShift - dPhiPart) / dPhiPart;
869 
870  fhPtHybrDetRecoil->Fill(ptbig,ptJet2);
871  fhPtHybrPartRecoil->Fill(ptbig,ptJet3);
872  fhPtHybrPartCorRecoil->Fill(ptcorr,ptJet3);
873  fhDPhiHybrPartCorRecoil->Fill(dPhiShift,dPhiPart);
874  fhResidualRecoil->Fill(residual);
875  fhPtResidualRecoil->Fill(ptJet3,residual);
876  fhDphiResidualRecoil->Fill(residualDphi);
877  fhDphiphiResidualRecoil->Fill(dPhiPart,residualDphi);
878 
879  if(fFillRecoilTree) {
880  fTreeVarsRecoil[0] = fCent;
881  fTreeVarsRecoil[1] = partback->Pt();
882  fTreeVarsRecoil[2] = ptbig;
883  fTreeVarsRecoil[3] = areabig;
884  fTreeVarsRecoil[4] = ptcorr;
885  fTreeVarsRecoil[5] = phibig;
886  fTreeVarsRecoil[6] = dPhiShift;
887  fTreeVarsRecoil[7] = ptJet3;
888  fTreeVarsRecoil[8] = phiJet3;
889  fTreeVarsRecoil[9] = dPhiPart;
891  fTreeEmbRecoil->Fill();
892  }
893  }
894  }
895  }
896 }
897 
899 
903  if(fDebug) Printf("particle container 0 entries = %i \t1 entries = %i\t 2 entries = %i",partCont0->GetNParticles(),partCont1->GetNParticles(),partCont2->GetNParticles());
907 
908  if(!jetCont || !jetContPart || !jetContTrue)
909  {
910  AliError(Form("jet container not found - check name %s",fJetContName.Data()));
911  TIter next(&fJetCollArray);
912  while ((jetCont = static_cast<AliJetContainer*>(next())))
913  AliError(Form("%s",jetCont->GetName()));
914  AliFatal("Exit...");
915  return;
916  }
917 
918  if(fDebug) {
919  Printf("n particle jets = %i",jetContPart->GetNJets());
920  Printf("n reco jets = %i",jetCont->GetNJets());
921  Printf("n PYTHIA jets = %i",jetContTrue->GetNJets());
922  }
923 
924  // Background
925  Double_t rho = 0;
926  if (jetCont->GetRhoParameter()) rho = jetCont->GetRhoVal();
927 
928  // PYTHIA event weight
929  // note - not used
930  // AliGenPythiaEventHeader *pyHeader = 0x0; //!<! Pythia header of the current external event
931  // AliAODEvent *ev = dynamic_cast<AliAODEvent*>(InputEvent());
932  // AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(ev->FindListObject(AliAODMCHeader::StdBranchName()));
933  // if (aodMCH) {
934  // for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
935  // pyHeader= dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
936  // if (pyHeader) break;
937  // }
938  // }
939  //
940  // Double_t pythiaCrossSection = 0;
941  // Double_t pythiaTrials = 0;
942  // Double_t pythiaWeight = 0;
943  // if (pyHeader)
944  // {
945  // if(fDebug) Printf("have pythia header - get weight");
946  // pythiaCrossSection = pyHeader->GetXsection();
947  // pythiaTrials = pyHeader->Trials();
948  // //fPythiaPtHard = fPythiaHeader->GetPtHard();
949  // pythiaWeight = pythiaCrossSection / pythiaTrials;
950  // }
951  // if(fDebug) Printf("pythia weight = %f",pythiaWeight);
952 
953  for(auto jet1 : jetCont->accepted()) { // loop over hybrid jets
954 
955  Double_t ptJet1 = jet1->Pt();
956  Double_t phiJet1 = jet1->Phi();
957  Double_t area = jet1->Area();
958  Double_t ptCorr = ptJet1-rho*area;
959  if(ptCorr<fMinEmbJetPt) continue;
960  auto jet2 = jet1->ClosestJet();
961  if(!jet2) {
962  //Printf("jet 2 cant be found");
963  continue;}
964  Double_t ptJet2 = jet2->Pt();
965  fhPtDet->Fill(ptJet2);
966  auto jet3 = jet2->ClosestJet();
967  if(!jet3) {
968  //Printf("jet3 can't be found");
969  continue;
970  }
971  fhPtDetMatchedToPart->Fill(ptJet2);
972  Double_t ptJet3 = jet3->Pt();
973  Double_t phiJet3 = jet3->Phi();
974 
975  if(phiJet3>2*TMath::Pi()) phiJet3 -= 2*TMath::Pi();
976  if(phiJet3<-2*TMath::Pi()) phiJet3 += 2*TMath::Pi();
977  if(phiJet3<-0.5*TMath::Pi()) phiJet3 += 2*TMath::Pi();
978  if(phiJet3>1.5*TMath::Pi()) phiJet3 -= 2*TMath::Pi();
979 
980  if(phiJet1>2*TMath::Pi()) phiJet1 -= 2*TMath::Pi();
981  if(phiJet1<-2*TMath::Pi()) phiJet1 += 2*TMath::Pi();
982  if(phiJet1<-0.5*TMath::Pi()) phiJet1 += 2*TMath::Pi();
983  if(phiJet1>1.5*TMath::Pi()) phiJet1 -= 2*TMath::Pi();
984 
985  if(fDebug) Printf("--- jet pt = jet hybrid pt = %f\t jet matched det pt = %f\t jet matched particle level pt = %f\t",ptJet1,ptJet2,ptJet3);
986 
987  fhPtDetPart->Fill(ptJet2,ptJet3);
988  Double_t fraction = jetCont->GetFractionSharedPt(jet1);
989  if(fraction < fMinFractionSharedPt) continue;
990 
991 
992  Double_t residual = (ptCorr - ptJet3) / ptJet3;
993  Double_t residualPhi = (phiJet1 - phiJet3) / phiJet3;
994 
995  fhPtHybrDet->Fill(ptJet1,ptJet2);
996  fhPtHybrPart->Fill(ptJet1,ptJet3);
997  fhPtHybrPartCor->Fill(ptCorr,ptJet3);
998  fhPhiHybrPartCor->Fill(phiJet1,phiJet3);
999 
1000  fhResidual->Fill(residual);
1001  fhPtResidual->Fill(ptJet3,residual);
1002  fhPhiResidual->Fill(residualPhi);
1003  fhPhiPhiResidual->Fill(phiJet3,residualPhi);
1004 
1007  fTreeVarsInclusive[1] = ptJet1;
1008  fTreeVarsInclusive[2] = area;
1009  fTreeVarsInclusive[3] = ptCorr;
1010  fTreeVarsInclusive[4] = phiJet1;
1011  fTreeVarsInclusive[5] = ptJet3;
1012  fTreeVarsInclusive[6] = phiJet3;
1014  fTreeEmbInclusive->Fill();
1015  }
1016  }
1017 // for(auto jettrue : jetContPart->accepted()) {
1018 // // jettrue
1019 // Double_t ptTrue = jettrue->Pt();
1020 // Double_t phiTrue = jettrue->Phi();
1021 // Double_t etaTrue= jettrue->Eta();
1022 // auto jetmatched = jettrue->ClosestJet();
1023 // if(!jetmatched) continue;
1024 // Double_t ptMatched = jetmatched->Pt();
1025 // Double_t phiMatched = jetmatched->Phi();
1026 // Double_t etaMatched = jetmatched->Eta();
1027 // Double_t residual = (ptMatched - ptTrue) / ptTrue;
1028 // fhPtHybrTrue->Fill(ptTrue,ptMatched);
1029 // fhResidual->Fill(residual);
1030 // fhPtResidual->Fill(ptTrue,residual);
1031 // }
1032 }
1033 
1034 
1035 
1041 {
1042  TString histname;
1043  TString groupname;
1044  AliJetContainer* jetCont = 0;
1045  TIter next(&fJetCollArray);
1046  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
1047  groupname = jetCont->GetName();
1048  UInt_t count = 0;
1049  for(auto jet : jetCont->accepted()) {
1050  if (!jet) continue;
1051  count++;
1052 
1053  histname = TString::Format("%s/histJetPt_%d", groupname.Data(), fCentBin);
1054  fHistManager.FillTH1(histname, jet->Pt());
1055 
1056  histname = TString::Format("%s/histJetArea_%d", groupname.Data(), fCentBin);
1057  fHistManager.FillTH1(histname, jet->Area());
1058 
1059  histname = TString::Format("%s/histJetPhi_%d", groupname.Data(), fCentBin);
1060  fHistManager.FillTH1(histname, jet->Phi());
1061 
1062  histname = TString::Format("%s/histJetEta_%d", groupname.Data(), fCentBin);
1063  fHistManager.FillTH1(histname, jet->Eta());
1064 
1065  if (jetCont->GetRhoParameter()) {
1066  histname = TString::Format("%s/histJetCorrPt_%d", groupname.Data(), fCentBin);
1067  fHistManager.FillTH1(histname, jet->Pt() - jetCont->GetRhoVal() * jet->Area());
1068  }
1069  }
1070  histname = TString::Format("%s/histNJets_%d", groupname.Data(), fCentBin);
1071  fHistManager.FillTH1(histname, count);
1072  }
1073 }
1074 
1080 {
1082 
1083  TString histname;
1084  TString groupname;
1085  UInt_t sumAcceptedTracks = 0;
1086  AliParticleContainer* partCont = 0;
1087  TIter next(&fParticleCollArray);
1088  while ((partCont = static_cast<AliParticleContainer*>(next()))) {
1089  groupname = partCont->GetName();
1090  UInt_t count = 0;
1091  for(auto part : partCont->accepted()) {
1092  if (!part) continue;
1093  count++;
1094 
1095  histname = TString::Format("%s/histTrackPt_%d", groupname.Data(), fCentBin);
1096  fHistManager.FillTH1(histname, part->Pt());
1097 
1098  histname = TString::Format("%s/histTrackPhi_%d", groupname.Data(), fCentBin);
1099  fHistManager.FillTH1(histname, part->Phi());
1100 
1101  histname = TString::Format("%s/histTrackEta_%d", groupname.Data(), fCentBin);
1102  fHistManager.FillTH1(histname, part->Eta());
1103 
1104  if (partCont->GetLoadedClass()->InheritsFrom("AliVTrack")) {
1105  const AliVTrack* track = static_cast<const AliVTrack*>(part);
1106 
1107  histname = TString::Format("%s/fHistDeltaEtaPt_%d", groupname.Data(), fCentBin);
1108  fHistManager.FillTH1(histname, track->Pt(), track->Eta() - track->GetTrackEtaOnEMCal());
1109 
1110  histname = TString::Format("%s/fHistDeltaPhiPt_%d", groupname.Data(), fCentBin);
1111  fHistManager.FillTH1(histname, track->Pt(), track->Phi() - track->GetTrackPhiOnEMCal());
1112 
1113  histname = TString::Format("%s/fHistDeltaPtvsPt_%d", groupname.Data(), fCentBin);
1114  fHistManager.FillTH1(histname, track->Pt(), track->Pt() - track->GetTrackPtOnEMCal());
1115 
1116  if (clusCont) {
1117  Int_t iCluster = track->GetEMCALcluster();
1118  if (iCluster >= 0) {
1119  AliVCluster* cluster = clusCont->GetAcceptCluster(iCluster);
1120  if (cluster) {
1121  histname = TString::Format("%s/fHistEoverPvsP_%d", groupname.Data(), fCentBin);
1122  fHistManager.FillTH2(histname, track->P(), cluster->GetNonLinCorrEnergy() / track->P());
1123  }
1124  }
1125  }
1126  }
1127  }
1128  sumAcceptedTracks += count;
1129 
1130  histname = TString::Format("%s/histNTracks_%d", groupname.Data(), fCentBin);
1131  fHistManager.FillTH1(histname, count);
1132  }
1133 
1134  histname = "fHistSumNTracks";
1135  fHistManager.FillTH1(histname, sumAcceptedTracks);
1136 }
1137 
1143 {
1144  TString histname;
1145  TString groupname;
1146  UInt_t sumAcceptedClusters = 0;
1147  AliClusterContainer* clusCont = 0;
1148  TIter next(&fClusterCollArray);
1149  while ((clusCont = static_cast<AliClusterContainer*>(next()))) {
1150  groupname = clusCont->GetName();
1151 
1152  for(auto cluster : clusCont->all()) {
1153  if (!cluster) continue;
1154 
1155  if (cluster->GetIsExotic()) {
1156  histname = TString::Format("%s/histClusterEnergyExotic_%d", groupname.Data(), fCentBin);
1157  fHistManager.FillTH1(histname, cluster->E());
1158  }
1159  }
1160 
1161  UInt_t count = 0;
1162  for(auto cluster : clusCont->accepted()) {
1163  if (!cluster) continue;
1164  count++;
1165 
1166  AliTLorentzVector nPart;
1167  cluster->GetMomentum(nPart, fVertex);
1168 
1169  histname = TString::Format("%s/histClusterEnergy_%d", groupname.Data(), fCentBin);
1170  fHistManager.FillTH1(histname, cluster->E());
1171 
1172  histname = TString::Format("%s/histClusterNonLinCorrEnergy_%d", groupname.Data(), fCentBin);
1173  fHistManager.FillTH1(histname, cluster->GetNonLinCorrEnergy());
1174 
1175  histname = TString::Format("%s/histClusterHadCorrEnergy_%d", groupname.Data(), fCentBin);
1176  fHistManager.FillTH1(histname, cluster->GetHadCorrEnergy());
1177 
1178  histname = TString::Format("%s/histClusterPhi_%d", groupname.Data(), fCentBin);
1179  fHistManager.FillTH1(histname, nPart.Phi_0_2pi());
1180 
1181  histname = TString::Format("%s/histClusterEta_%d", groupname.Data(), fCentBin);
1182  fHistManager.FillTH1(histname, nPart.Eta());
1183  }
1184  sumAcceptedClusters += count;
1185 
1186  histname = TString::Format("%s/histNClusters_%d", groupname.Data(), fCentBin);
1187  fHistManager.FillTH1(histname, count);
1188  }
1189 
1190  histname = "fHistSumNClusters";
1191  fHistManager.FillTH1(histname, sumAcceptedClusters);
1192 }
1193 
1199 {
1200  if (!fCaloCells) return;
1201 
1202  TString histname;
1203 
1204  const Short_t ncells = fCaloCells->GetNumberOfCells();
1205 
1206  histname = TString::Format("%s/histNCells_%d", fCaloCellsName.Data(), fCentBin);
1207  fHistManager.FillTH1(histname, ncells);
1208 
1209  histname = TString::Format("%s/histCellEnergy_%d", fCaloCellsName.Data(), fCentBin);
1210  for (Short_t pos = 0; pos < ncells; pos++) {
1211  Double_t amp = fCaloCells->GetAmplitude(pos);
1212 
1213  fHistManager.FillTH1(histname, amp);
1214  }
1215 }
1216 
1222 {
1224 }
1225 
1234 {
1235  return kTRUE;
1236 }
1237 
1242 {
1243 }
1244 
1246 
1247  Int_t index=-1;
1248  Int_t triggers[100];
1249 
1250  for(Int_t cr=0;cr<100;cr++) triggers[cr]=-1;
1251 
1252  Int_t im=0;
1253 
1254  TString groupname = "";
1255  AliParticleContainer* partCont = 0x0;
1257  else partCont = GetParticleContainer(0);
1258  groupname = partCont->GetName();
1259  UInt_t iCount = 0;
1260  for(auto part : partCont->accepted()) {
1261  if (!part) continue;
1262  list->Add(part);
1263  iCount++;
1264  if(part->Pt()>=minT && part->Pt()<maxT){
1265  triggers[im]=iCount-1;
1266  im=im+1;
1267  }
1268  }
1269  number=im;
1270  Int_t rd=0;
1271  if(im>0) rd=fRandom->Integer(im);
1272  index=triggers[rd];
1273  return index;
1274 }
1275 
1277 
1278  // get relative DeltaPhi from -pi < DeltaPhi < pi
1279 
1280  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1281  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1282  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1283  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1284  double dphi = mphi-vphi;
1285  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1286  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1287  return dphi;//dphi in [-Pi, Pi]
1288 }
1289 
1291 {
1292  Int_t phibin=-1;
1293  if(!(TMath::Abs(phi)<=2*TMath::Pi())) return -1;
1294  Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1295  phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1296  //if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1297  return phibin;
1298 }
Float_t fJetHadronDeltaPhi
max angle from pi (set <0 for no selection)
THashList * CreateHistoGroup(const char *groupname)
Create a new group of histograms within a parent group.
TObjArray fClusterCollArray
cluster collection array
Float_t fTTLowRef
minimum reference trigger track pt
Double_t GetRhoVal() const
const TString & GetRhoName() const
double Double_t
Definition: External.C:58
Float_t fFrac
fraction of events that are used to fill signal recoil jet population
Bool_t fFillRecoilTHnSparse
switch to fill recoil THnSparse for main analysis
Definition: External.C:236
Int_t SelectTrigger(TList *list, Double_t minT, Double_t maxT, Int_t &number)
AliJetContainer * GetJetContainer(Int_t i=0) const
Float_t fTTUpRef
maximum reference trigger track pt
Int_t GetNParticles() const
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
Declaration of class AliTLorentzVector.
Bool_t fFillInclusiveTree
switch to fill embedding tree with inclusive jet info
Double_t fMinBinPt
min pt in histograms
Double_t fEPV0
!event plane V0
Int_t fCentBin
!event centrality bin
Float_t fMinEmbJetPt
min corrected jet pt to use in embedding
TString fJetContName
Base level jet container name.
Float_t fMinFractionSharedPt
min fraction of pt between hybrid / detector jets
Container for particles within the EMCAL framework.
TString fJetContTrueName
True pp (detector) level jet container name.
TString kData
Declare data MC or deltaAOD.
TObjArray fParticleCollArray
particle/track collection array
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
const AliClusterIterableContainer all() const
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
Create a new TH2 within the container.
TObject * FindObject(const char *name) const
Find an object inside the container.
int Int_t
Definition: External.C:63
Float_t fTTUpSig
maximum signal trigger track pt
Definition: External.C:204
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
float Float_t
Definition: External.C:68
Double_t Phi_0_2pi() const
Float_t fTTLowSig
minimum signal trigger track pt
Int_t GetNJets() const
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
BeamType fForceBeamType
forced beam type
Int_t fNcentBins
how many centrality bins
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
Definition: External.C:212
AliVCluster * GetAcceptCluster(Int_t i) const
TString fJetContPartName
Particle(MC) level jet container name.
const AliClusterIterableContainer accepted() const
Double_t fCent
!event centrality
Float_t fCentMax
maximum centrality
TString fCaloCellsName
name of calo cell collection
Bool_t fFillJetHistograms
switch to fill jet histograms
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
TObjArray fJetCollArray
jet collection array
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Double_t RelativePhi(Double_t mphi, Double_t vphi)
short Short_t
Definition: External.C:23
AliVCaloCells * fCaloCells
!cells
AliRhoParameter * GetRhoParameter()
AliEmcalList * fOutput
!output list
Double_t fMaxBinPt
max pt in histograms
Int_t fPtHardBin
pt hard bin if running embedding
THistManager fHistManager
Histogram manager.
Double_t fVertex[3]
!event vertex
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
Double_t GetFractionSharedPt(const AliEmcalJet *jet, AliParticleContainer *cont2=0x0) const
Declaration of class AliEmcalPythiaInfo.
const AliParticleIterableContainer accepted() const
const char Option_t
Definition: External.C:48
void UserCreateOutputObjects()
Main initialization function on the worker.
const AliJetIterableContainer accepted() const
bool Bool_t
Definition: External.C:53
Bool_t fFillRecoilTree
switch to fill embedding tree with recoil jet info
Container structure for EMCAL clusters.
Float_t fCentMin
minimum centrality
Container for jet within the EMCAL jet framework.
Definition: External.C:196
Bool_t fFillTrackHistograms
switch to fill track histograms
Int_t fNbins
no. of pt bins
Int_t fRejectionFactorInclusiveJets
factor to reject inclusive jets, to reduce size of ttree