AliPhysics  e6c8d43 (e6c8d43)
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  fJetEtaMin(-.5),
68  fJetEtaMax(.5),
69  fJetHadronDeltaPhi(0.6),
70  fMinFractionSharedPt(0.5),
71  fJetContName(""),
72  fJetContTrueName(""),
73  fJetContPartName(""),
74  fFillTrackHistograms(kTRUE),
75  fFillJetHistograms(kTRUE),
76  fFillInclusiveTree(kFALSE),
77  fFillRecoilTree(kFALSE),
78  fPtHardBin(0.),
79  fRandom(0),
80  fHistEvtSelection(0x0),
81  fHJetSpec(0x0),
82  fh1TrigRef(0x0),
83  fh1TrigSig(0x0),
84  fh2Ntriggers(0x0),
85  fh2RPJetsC10(0x0),
86  fh2RPJetsC20(0x0),
87  fh2RPTC10(0x0),
88  fh2RPTC20(0x0),
89  fhDphiPtSig(0x0),
90  fhDphiPtRef(0x0),
91  fhPtDetPart(0x0),
92  fhPtHybrDet(0x0),
93  fhPtHybrPart(0x0),
94  fhPtHybrPartCor(0x0),
95  fhPhiHybrPartCor(0x0),
96  fhPtDet(0x0),
97  fhPtDetMatchedToPart(0x0),
98  fhResidual(0x0),
99  fhPtResidual(0x0),
100  fhPhiResidual(0x0),
101  fhPhiPhiResidual(0x0),
102  fhPtDetPartRecoil(0x0),
103  fhPtHybrDetRecoil(0x0),
104  fhPtHybrPartRecoil(0x0),
105  fhPtHybrPartCorRecoil(0x0),
106  fhPtDetRecoil(0x0),
107  fhPtDetMatchedToPartRecoil(0x0),
108  fhResidualRecoil(0x0),
109  fhPtResidualRecoil(0x0),
110  fhDphiResidualRecoil(0x0),
111  fhDphiphiResidualRecoil(0x0),
112  fhTTPtDetMatchedToPart(0x0),
113  fhTTPhiDetMatchedToPart(0x0),
114  fhDPhiHybrPartCorRecoil(0x0),
115  fTreeEmbInclusive(0x0),
116  fTreeEmbRecoil(0x0)
117 {
119 
120  DefineOutput(1, TList::Class());
121  //if(fJetShapeType == AliAnalysisTaskJetCoreEmcal::kDetEmbPart && fFillInclusiveTree)
122  DefineOutput(2, TTree::Class());
123 // if(fJetShapeType == AliAnalysisTaskJetCoreEmcal::kDetEmbPart && fFillRecoilTree)
124  DefineOutput(3, TTree::Class());
125 }
126 
133  AliAnalysisTaskEmcalJet(name, kTRUE),
134  fHistManager(name),
136  fCentMin(0.),
137  fCentMax(100.),
138  fTTLowRef(8),
139  fTTUpRef(9.),
140  fTTLowSig(20.),
141  fTTUpSig(50.),
142  fNRPBins(50),
143  fFrac(0.8),
144  fJetEtaMin(-.5),
145  fJetEtaMax(.5),
146  fJetHadronDeltaPhi(0.6),
148  fJetContName(""),
149  fJetContTrueName(""),
150  fJetContPartName(""),
151  fFillTrackHistograms(kTRUE),
152  fFillJetHistograms(kTRUE),
153  fFillInclusiveTree(kFALSE),
154  fFillRecoilTree(kFALSE),
155  fPtHardBin(0.),
156  fRandom(0),
157  fHistEvtSelection(0x0),
158  fHJetSpec(0x0),
159  fh1TrigRef(0x0),
160  fh1TrigSig(0x0),
161  fh2Ntriggers(0x0),
162  fh2RPJetsC10(0x0),
163  fh2RPJetsC20(0x0),
164  fh2RPTC10(0x0),
165  fh2RPTC20(0x0),
166  fhDphiPtSig(0x0),
167  fhDphiPtRef(0x0),
168  fhPtDetPart(0x0),
169  fhPtHybrDet(0x0),
170  fhPtHybrPart(0x0),
171  fhPtHybrPartCor(0x0),
172  fhPhiHybrPartCor(0x0),
173  fhPtDet(0x0),
175  fhResidual(0x0),
176  fhPtResidual(0x0),
177  fhPhiResidual(0x0),
178  fhPhiPhiResidual(0x0),
179  fhPtDetPartRecoil(0x0),
180  fhPtHybrDetRecoil(0x0),
181  fhPtHybrPartRecoil(0x0),
183  fhPtDetRecoil(0x0),
185  fhResidualRecoil(0x0),
186  fhPtResidualRecoil(0x0),
192  fTreeEmbInclusive(0x0),
193  fTreeEmbRecoil(0x0)
194 {
196 
197  DefineOutput(1, TList::Class());
198  DefineOutput(2, TTree::Class());
199  DefineOutput(3, TTree::Class());
200 }
201 
206 {
207 }
208 
214 {
216 
219  //AllocateClusterHistograms();
220  //AllocateCellHistograms();
222 
223  TIter next(fHistManager.GetListOfHistograms());
224  TObject* obj = 0;
225  while ((obj = next())) {
226  fOutput->Add(obj);
227  }
228 
229  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
230  if(fJetShapeType == AliAnalysisTaskJetCoreEmcal::kDetEmbPart && fFillInclusiveTree) PostData(2, fTreeEmbInclusive); // Post data for ALL output slots > 0 here.
231  if(fJetShapeType == AliAnalysisTaskJetCoreEmcal::kDetEmbPart && fFillRecoilTree) PostData(3, fTreeEmbRecoil); // Post data for ALL output slots > 0 here.
232 }
233 
234 /*
235  * This function allocates the histograms for basic EMCal cluster QA.
236  * A set of histograms (energy, eta, phi, number of cluster) is allocated
237  * per each cluster container and per each centrality bin.
238  */
240 {
241  TString histname;
242  TString histtitle;
243  TString groupname;
244  AliClusterContainer* clusCont = 0;
245  TIter next(&fClusterCollArray);
246  while ((clusCont = static_cast<AliClusterContainer*>(next()))) {
247  groupname = clusCont->GetName();
248  // Protect against creating the histograms twice
249  if (fHistManager.FindObject(groupname)) {
250  AliWarning(TString::Format("%s: Found groupname %s in hist manager. The cluster containers will be filled into the same histograms.", GetName(), groupname.Data()));
251  continue;
252  }
253  fHistManager.CreateHistoGroup(groupname);
254  for (Int_t cent = 0; cent < fNcentBins; cent++) {
255  histname = TString::Format("%s/histClusterEnergy_%d", groupname.Data(), cent);
256  histtitle = TString::Format("%s;#it{E}_{cluster} (GeV);counts", histname.Data());
257  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
258 
259  histname = TString::Format("%s/histClusterEnergyExotic_%d", groupname.Data(), cent);
260  histtitle = TString::Format("%s;#it{E}_{cluster}^{exotic} (GeV);counts", histname.Data());
261  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
262 
263  histname = TString::Format("%s/histClusterNonLinCorrEnergy_%d", groupname.Data(), cent);
264  histtitle = TString::Format("%s;#it{E}_{cluster}^{non-lin.corr.} (GeV);counts", histname.Data());
265  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
266 
267  histname = TString::Format("%s/histClusterHadCorrEnergy_%d", groupname.Data(), cent);
268  histtitle = TString::Format("%s;#it{E}_{cluster}^{had.corr.} (GeV);counts", histname.Data());
269  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
270 
271  histname = TString::Format("%s/histClusterPhi_%d", groupname.Data(), cent);
272  histtitle = TString::Format("%s;#it{#phi}_{custer};counts", histname.Data());
273  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
274 
275  histname = TString::Format("%s/histClusterEta_%d", groupname.Data(), cent);
276  histtitle = TString::Format("%s;#it{#eta}_{custer};counts", histname.Data());
277  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
278 
279  histname = TString::Format("%s/histNClusters_%d", groupname.Data(), cent);
280  histtitle = TString::Format("%s;number of clusters;events", histname.Data());
281  if (fForceBeamType != kpp) {
282  fHistManager.CreateTH1(histname, histtitle, 500, 0, 3000);
283  }
284  else {
285  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
286  }
287  }
288  }
289 
290  histname = "fHistSumNClusters";
291  histtitle = TString::Format("%s;Sum of n clusters;events", histname.Data());
292  if (fForceBeamType != kpp) {
293  fHistManager.CreateTH1(histname, histtitle, 500, 0, 3000);
294  }
295  else {
296  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
297  }
298 }
299 
300 /*
301  * This function allocates the histograms for basic EMCal QA.
302  * One 2D histogram with the cell energy spectra and the number of cells
303  * per event is allocated per each centrality bin.
304  */
306 {
307  TString histname;
308  TString histtitle;
309  TString groupname(fCaloCellsName);
310 
311  fHistManager.CreateHistoGroup(groupname);
312  for (Int_t cent = 0; cent < fNcentBins; cent++) {
313  histname = TString::Format("%s/histCellEnergy_%d", groupname.Data(), cent);
314  histtitle = TString::Format("%s;#it{E}_{cell} (GeV);counts", histname.Data());
315  fHistManager.CreateTH1(histname, histtitle, 300, 0, 150);
316 
317  histname = TString::Format("%s/histNCells_%d", groupname.Data(), cent);
318  histtitle = TString::Format("%s;number of cells;events", histname.Data());
319  if (fForceBeamType != kpp) {
320  fHistManager.CreateTH1(histname, histtitle, 500, 0, 6000);
321  }
322  else {
323  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
324  }
325  }
326 }
327 
328 /*
329  * This function allocates the histograms for basic tracking QA.
330  * A set of histograms (pT, eta, phi, difference between kinematic properties
331  * at the vertex and at the EMCal surface, number of tracks) is allocated
332  * per each particle container and per each centrality bin.
333  */
335 {
336  TString histname;
337  TString histtitle;
338  TString groupname;
339  AliParticleContainer* partCont = 0;
340  TIter next(&fParticleCollArray);
341  while ((partCont = static_cast<AliParticleContainer*>(next()))) {
342  groupname = partCont->GetName();
343 
344  // Protect against creating the histograms twice
345  if (fHistManager.FindObject(groupname)) {
346  AliWarning(TString::Format("%s: Found groupname %s in hist manager. The track containers will be filled into the same histograms.", GetName(), groupname.Data()));
347  continue;
348  }
349  fHistManager.CreateHistoGroup(groupname);
350  for (Int_t cent = 0; cent < fNcentBins; cent++) {
351  histname = TString::Format("%s/histTrackPt_%d", groupname.Data(), cent);
352  histtitle = TString::Format("%s;#it{p}_{T,track} (GeV/#it{c});counts", histname.Data());
353  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
354 
355  histname = TString::Format("%s/histTrackPhi_%d", groupname.Data(), cent);
356  histtitle = TString::Format("%s;#it{#phi}_{track};counts", histname.Data());
357  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
358 
359  histname = TString::Format("%s/histTrackEta_%d", groupname.Data(), cent);
360  histtitle = TString::Format("%s;#it{#eta}_{track};counts", histname.Data());
361  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
362 
363  if (TClass(partCont->GetClassName()).InheritsFrom("AliVTrack")) {
364  histname = TString::Format("%s/fHistDeltaEtaPt_%d", groupname.Data(), cent);
365  histtitle = TString::Format("%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{#eta}_{track}^{vertex} - #it{#eta}_{track}^{EMCal};counts", histname.Data());
366  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, 50, -0.5, 0.5);
367 
368  histname = TString::Format("%s/fHistDeltaPhiPt_%d", groupname.Data(), cent);
369  histtitle = TString::Format("%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{#phi}_{track}^{vertex} - #it{#phi}_{track}^{EMCal};counts", histname.Data());
370  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, 200, -2, 2);
371 
372  histname = TString::Format("%s/fHistDeltaPtvsPt_%d", groupname.Data(), cent);
373  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());
374  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, fNbins / 2, -fMaxBinPt/2, fMaxBinPt/2);
375 
376  histname = TString::Format("%s/fHistEoverPvsP_%d", groupname.Data(), cent);
377  histtitle = TString::Format("%s;#it{P}_{track} (GeV/#it{c});#it{E}_{cluster} / #it{P}_{track} #it{c};counts", histname.Data());
378  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, fNbins / 2, 0, 4);
379  }
380 
381  histname = TString::Format("%s/histNTracks_%d", groupname.Data(), cent);
382  histtitle = TString::Format("%s;number of tracks;events", histname.Data());
383  if (fForceBeamType != kpp) {
384  fHistManager.CreateTH1(histname, histtitle, 500, 0, 5000);
385  }
386  else {
387  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
388  }
389  }
390  }
391 
392  histname = "fHistSumNTracks";
393  histtitle = TString::Format("%s;Sum of n tracks;events", histname.Data());
394  if (fForceBeamType != kpp) {
395  fHistManager.CreateTH1(histname, histtitle, 500, 0, 5000);
396  }
397  else {
398  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
399  }
400 }
401 
402 /*
403  * This function allocates the histograms for basic jet QA.
404  * A set of histograms (pT, eta, phi, area, number of jets, corrected pT) is allocated
405  * per each jet container and per each centrality bin.
406  */
408 {
409  TString histname;
410  TString histtitle;
411  TString groupname;
412  AliJetContainer* jetCont = 0;
413  TIter next(&fJetCollArray);
414  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
415  groupname = jetCont->GetName();
416  // Protect against creating the histograms twice
417  if (fHistManager.FindObject(groupname)) {
418  AliWarning(TString::Format("%s: Found groupname %s in hist manager. The jet containers will be filled into the same histograms.", GetName(), groupname.Data()));
419  continue;
420  }
421  fHistManager.CreateHistoGroup(groupname);
422  for (Int_t cent = 0; cent < fNcentBins; cent++) {
423  histname = TString::Format("%s/histJetPt_%d", groupname.Data(), cent);
424  histtitle = TString::Format("%s;#it{p}_{T,jet} (GeV/#it{c});counts", histname.Data());
425  fHistManager.CreateTH1(histname, histtitle, fNbins, fMinBinPt, fMaxBinPt);
426 
427  histname = TString::Format("%s/histJetArea_%d", groupname.Data(), cent);
428  histtitle = TString::Format("%s;#it{A}_{jet};counts", histname.Data());
429  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, 3);
430 
431  histname = TString::Format("%s/histJetPhi_%d", groupname.Data(), cent);
432  histtitle = TString::Format("%s;#it{#phi}_{jet};counts", histname.Data());
433  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
434 
435  histname = TString::Format("%s/histJetEta_%d", groupname.Data(), cent);
436  histtitle = TString::Format("%s;#it{#eta}_{jet};counts", histname.Data());
437  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
438 
439  histname = TString::Format("%s/histNJets_%d", groupname.Data(), cent);
440  histtitle = TString::Format("%s;number of jets;events", histname.Data());
441  if (fForceBeamType != kpp) {
442  fHistManager.CreateTH1(histname, histtitle, 500, 0, 500);
443  }
444  else {
445  fHistManager.CreateTH1(histname, histtitle, 100, 0, 100);
446  }
447 
448  if (!jetCont->GetRhoName().IsNull()) {
449  histname = TString::Format("%s/histJetCorrPt_%d", groupname.Data(), cent);
450  histtitle = TString::Format("%s;#it{p}_{T,jet}^{corr} (GeV/#it{c});counts", histname.Data());
451  fHistManager.CreateTH1(histname, histtitle, fNbins, -fMaxBinPt / 2, fMaxBinPt / 2);
452  }
453  }
454  }
455 }
456 
458 {
459 
460  Bool_t oldStatus = TH1::AddDirectoryStatus();
461  TH1::AddDirectory(kFALSE);
462 
463  // set seed
464  fRandom = new TRandom3(0);
465 
466  fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
467  fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
468  fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
469  fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
470  fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
471  fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
472  fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
473 
474  fh1TrigRef=new TH1D("Trig Ref","",10,0.,10);
475  fh1TrigSig=new TH1D("Trig Sig","",10,0.,10);
476  fh2Ntriggers=new TH2F("# of triggers","",100,0.,100.,50,0.,50.);
477  fh2RPJetsC10=new TH2F("RPJetC10","",35,0.,3.5,100,0.,100.);
478  fh2RPJetsC20=new TH2F("RPJetC20","",35,0.,3.5,100,0.,100.);
479  fh2RPTC10=new TH2F("RPTriggerC10","",35,0.,3.5,50,0.,50.);
480  fh2RPTC20=new TH2F("RPTriggerC20","",35,0.,3.5,50,0.,50.);
481 
483 
484  fOutput->Add(fh1TrigRef);
485  fOutput->Add(fh1TrigSig);
486  fOutput->Add(fh2Ntriggers);
487  fOutput->Add(fh2RPJetsC10);
488  fOutput->Add(fh2RPJetsC20);
489  fOutput->Add(fh2RPTC10);
490  fOutput->Add(fh2RPTC20);
491 
492  const Int_t dimSpec = 6;
493  const Int_t nBinsSpec[dimSpec] = {100,100, 280, 50,200, fNRPBins};
494  const Double_t lowBinSpec[dimSpec] = {0,0,-80, 0,-0.5*TMath::Pi(), 0};
495  const Double_t hiBinSpec[dimSpec] = {100,1, 200, 50,1.5*TMath::Pi(), static_cast<Double_t>(fNRPBins)};
496  fHJetSpec = new THnSparseF("fHJetSpec","Recoil jet spectrum",dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
497 
498  // comment out since I want finer binning in jet area, to make it easier
499  // to change selection on jet area (Leticia used 0.8*R^2*Pi whereas 0.6 is used
500  // for inclusive jets)
501  //change binning in jet area
502 // Double_t *xPt6 = new Double_t[7];
503 // xPt6[0] = 0.;
504 // xPt6[1]=0.07;
505 // xPt6[2]=0.2;
506 // xPt6[3]=0.4;
507 // xPt6[4]=0.6;
508 // xPt6[5]=0.8;
509 // xPt6[6]=1;
510 // fHJetSpec->SetBinEdges(1,xPt6);
511 // delete [] xPt6;
512 
513  fOutput->Add(fHJetSpec);
514 
515  // azimuthal correlation
516 
517  fhDphiPtSig = new TH2F("hDphiPtS","recoil #Delta #phi vs jet pT signal",100,-2,5,250,-50,200);
518  fhDphiPtSig->GetXaxis()->SetTitle("#Delta #phi");
519  fhDphiPtSig->GetYaxis()->SetTitle("p^{reco,ch}_{T,jet} (GeV/c)");
520  fhDphiPtRef = new TH2F("hDphiPtR","recoil #Delta #phi vs jet pT reference",100,-2,5,250,-50,200);
521  fhDphiPtRef->GetXaxis()->SetTitle("#Delta #phi");
522  fhDphiPtRef->GetYaxis()->SetTitle("p^{reco,ch}_{T,jet} (GeV/c)");
523 
524  fOutput->Add(fhDphiPtRef);
525  fOutput->Add(fhDphiPtSig);
526 
527 
528  fhPtHybrDet= new TH2F("hPtHybrDet","pT response Pb-Pb+PYTHIA vs PYTHIA",200,0,200,200,0,200);
529  fhPtHybrDet->GetXaxis()->SetTitle("p^{Pb-Pb+PYTHIA,ch}_{T} (GeV/c)");
530  fhPtHybrDet->GetYaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
531  fhPtDetPart= new TH2F("hPtDetPart","pT response PYTHIA vs part",200,0,200,200,0,200);
532  fhPtDetPart->GetXaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
533  fhPtDetPart->GetYaxis()->SetTitle("p^{true}_{T} (GeV/c)");
534  fhPtHybrPart = new TH2F("hPtHybrPart",Form("pT response Pb-Pb+PYTHIA vs part, min shared pT > %.0f",fMinFractionSharedPt*100),200,0,200,200,0,200);
535  fhPtHybrPart->GetXaxis()->SetTitle("p^{Pb-Pb+PYTHIA,ch}_{T} (GeV/c)");
536  fhPtHybrPart->GetYaxis()->SetTitle("p^{true}_{T} (GeV/c)");
537  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);
538  fhPtHybrPartCor->GetXaxis()->SetTitle("p^{Pb-Pb+PYTHIA,ch}_{T} (GeV/c)");
539  fhPtHybrPartCor->GetYaxis()->SetTitle("p^{true}_{T} (GeV/c)");
540 
541  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());
542  fhPhiHybrPartCor->GetXaxis()->SetTitle("#phi^{Pb-Pb+PYTHIA,ch}");
543  fhPhiHybrPartCor->GetYaxis()->SetTitle("#phi^{true}");
544  fhResidual = new TH1F("hResidual","residual",50,-1,1);
545  fhResidual->GetXaxis()->SetTitle("p^{reco}_{T} - p^{true}_{T} / p^{true}_{T}");
546  fhPtResidual= new TH2F("hPtResidual","pT vs residual",200,0,200,50,-1,1);
547  fhPtResidual->GetXaxis()->SetTitle("p^{true}_{T} (GeV/c)");
548  fhPtResidual->GetYaxis()->SetTitle("p^{reco}_{T} - p^{true}_{T} / p^{true}_{T}");
549  fhPhiResidual = new TH1F("hPhiResidual","residual phi",50,-1,1);
550  fhPhiResidual->GetXaxis()->SetTitle("#phi^{reco} - #phi^{true} / #phi^{true}");
551  fhPhiPhiResidual= new TH2F("hPhiPhiResidual","pT vs residual",600,-3,3,200,-0.5*TMath::Pi(),1.5*TMath::Pi());
552  fhPhiPhiResidual->GetXaxis()->SetTitle("#phi^{true}");
553  fhPhiPhiResidual->GetYaxis()->SetTitle("#phi^{reco} - #phi^{true} / #phi^{true}");
554 
555  fhPtDet= new TH1F("hPtDet","pT detector level",200,0,200);
556  fhPtDet->GetXaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
557  fhPtDetMatchedToPart = new TH1F("hPtDetMatchedToPart","pT detector level matched to particle level jet",200,0,200);
558  fhPtDetMatchedToPart->GetXaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
559 
560  fhPtHybrDetRecoil= new TH2F("hPtHybrDetRecoil","pT response Pb-Pb+PYTHIA vs PYTHIA",200,0,200,200,0,200);
561  fhPtHybrDetRecoil->GetXaxis()->SetTitle("p^{Pb-Pb+PYTHIA,ch}_{T} (GeV/c)");
562  fhPtHybrDetRecoil->GetYaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
563  fhPtDetPartRecoil= new TH2F("hPtDetPartRecoil","pT response PYTHIA vs part",200,0,200,200,0,200);
564  fhPtDetPartRecoil->GetXaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
565  fhPtDetPartRecoil->GetYaxis()->SetTitle("p^{true}_{T} (GeV/c)");
566  fhPtHybrPartRecoil = new TH2F("hPtHybrPartRecoil",Form("pT response Pb-Pb+PYTHIA vs part, min shared pT > %.0f",fMinFractionSharedPt*100),200,0,200,200,0,200);
567  fhPtHybrPartRecoil->GetXaxis()->SetTitle("p^{Pb-Pb+PYTHIA,ch}_{T} (GeV/c)");
568  fhPtHybrPartRecoil->GetYaxis()->SetTitle("p^{true}_{T} (GeV/c)");
569  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);
570  fhPtHybrPartCorRecoil->GetXaxis()->SetTitle("p^{Pb-Pb+PYTHIA,ch}_{T} (GeV/c)");
571  fhPtHybrPartCorRecoil->GetYaxis()->SetTitle("p^{true}_{T} (GeV/c)");
572 
573  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());
574  fhDPhiHybrPartCorRecoil->GetXaxis()->SetTitle("#Delta#phi^{Pb-Pb+PYTHIA,ch}");
575  fhDPhiHybrPartCorRecoil->GetYaxis()->SetTitle("#Delta#phi^{true}");
576  fhResidualRecoil = new TH1F("hResidualRecoil","residual",50,-1,1);
577  fhResidualRecoil->GetXaxis()->SetTitle("p^{reco}_{T} - p^{true}_{T} / p^{true}_{T}");
578  fhPtResidualRecoil= new TH2F("hPtResidualRecoil","pT vs residual",200,0,200,50,-1,1);
579  fhPtResidualRecoil->GetXaxis()->SetTitle("p^{true}_{T} (GeV/c)");
580  fhPtResidualRecoil->GetYaxis()->SetTitle("p^{reco}_{T} - p^{true}_{T} / p^{true}_{T}");
581 
582  fhDphiResidualRecoil = new TH1F("hDphiResidualRecoil","residual",50,-1,1);
583  fhDphiResidualRecoil->GetXaxis()->SetTitle("#Delta#phi^{reco} - #Delta#phi^{true} / #Delta#phi^{true}");
584  fhDphiphiResidualRecoil= new TH2F("hDphiphiResidualRecoil","#Delta#phi vs residual in phi",400,-3,3,200,-0.5*TMath::Pi(),1.5*TMath::Pi());
585  fhDphiphiResidualRecoil->GetXaxis()->SetTitle("#Delta#phi^{true}");
586  fhDphiphiResidualRecoil->GetYaxis()->SetTitle("#Delta#phi^{reco} - #Delta#phi^{true} / #Delta#phi^{true}");
587  fhPtDetRecoil= new TH1F("hPtDetRecoil","pT detector level",200,0,200);
588  fhPtDetRecoil->GetXaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
589  fhPtDetMatchedToPartRecoil = new TH1F("hPtDetMatchedToPartRecoil","pT detector level matched to particle level jet",200,0,200);
590  fhPtDetMatchedToPartRecoil->GetXaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
591 
592  fhTTPtDetMatchedToPart = new TH2F("hTTPtDetMatchedToPart","trigger track pT response reco vs partice",140,0,70,140,0,70);
593  fhTTPtDetMatchedToPart->GetXaxis()->SetTitle("p^{TT,reco}_{T} (GeV/c)");
594  fhTTPtDetMatchedToPart->GetYaxis()->SetTitle("p^{TT,part}_{T} (GeV/c)");
595  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());
596  fhTTPhiDetMatchedToPart->GetXaxis()->SetTitle("p^{TT,reco}_{T} (GeV/c)");
597  fhTTPhiDetMatchedToPart->GetYaxis()->SetTitle("p^{TT,part}_{T} (GeV/c)");
598 
599  fOutput->Add(fhPtHybrDet);
600  fOutput->Add(fhPtDetPart);
601  fOutput->Add(fhPtHybrPart);
602  fOutput->Add(fhPtHybrPartCor);
604  fOutput->Add(fhPtDet);
606  fOutput->Add(fhResidual);
607  fOutput->Add(fhPtResidual);
608  fOutput->Add(fhPhiResidual);
610 
616  fOutput->Add(fhPtDetRecoil);
622 
625 
626  TString varNamesInclusive[8]={"centrality","ptRawRec","areaRec","ptCorrRec","phiRec","ptPart","phiPart","binPtHard"};
627  TString varNamesRecoil[11]={"centrality","ptTT","ptRawRec","areaRec","ptCorrRec","phiRec","DPhiRec","ptPart","phiPart","DPhiPart","binPtHard"};
629  const char* nameEmbInclusive = GetOutputSlot(2)->GetContainer()->GetName();
630  fTreeEmbInclusive = new TTree(nameEmbInclusive, nameEmbInclusive);
631  for(Int_t ivar=0; ivar < 8; ivar++){
632  fTreeEmbInclusive->Branch(varNamesInclusive[ivar].Data(), &fTreeVarsInclusive[ivar], Form("%s/F", varNamesInclusive[ivar].Data()));
633  }
634  }
635 
637  const char* nameEmbRecoil= GetOutputSlot(3)->GetContainer()->GetName();
638  fTreeEmbRecoil = new TTree(nameEmbRecoil, nameEmbRecoil);
639  for(Int_t ivar=0; ivar < 11; ivar++){
640  fTreeEmbRecoil->Branch(varNamesRecoil[ivar].Data(), &fTreeVarsRecoil[ivar], Form("%s/F", varNamesRecoil[ivar].Data()));
641  }
642  }
643 
644  // =========== Switch on Sumw2 for all histos ===========
645  for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
646  TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
647  if (h1){
648  h1->Sumw2();
649  continue;
650  }
651  THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
652  if (hn){
653  hn->Sumw2();
654  }
655  }
656 
657  // add QA plots from fEventCuts
658  fEventCuts.AddQAplotsToList(fOutput);
659 
660  TH1::AddDirectory(oldStatus);
661 
662  PostData(1, fOutput);
665 }
666 
674 {
675 
676  fHistEvtSelection->Fill(1); // number of events before event selection
677  AliVEvent *ev = InputEvent();
678  if (!fEventCuts.AcceptEvent(ev)) {
679  fHistEvtSelection->Fill(2);
680  return kTRUE;
681  }
682 
685  //DoClusterLoop();
686  //DoCellLoop();
687  DoJetCoreLoop();
689 
690  return kTRUE;
691 }
692 
694 {
695  // Do jet core analysis and fill histograms.
696 
698  AliJetContainer *jetContPart = 0x0;
699  AliJetContainer *jetContTrue = 0x0;
701  jetContTrue = GetJetContainer(fJetContTrueName);
702  jetContPart = GetJetContainer(fJetContPartName);
703  }
704 
705  if(!jetCont ||
707  AliError(Form("jet container not found - check name %s",fJetContName.Data()));
708  TIter next(&fJetCollArray);
709  while ((jetCont = static_cast<AliJetContainer*>(next())))
710  AliError(Form("%s",jetCont->GetName()));
711  AliFatal("Exit...");
712  return;
713  }
714 
715  // centrality selection
716  if(fDebug) Printf("centrality: %f\n", fCent);
718  ((fCent>fCentMax) || (fCent<fCentMin))) {
719  fHistEvtSelection->Fill(4);
720  return;
721  }
722  fHistEvtSelection->Fill(0);
723 
724  // Background
725  Double_t rho = 0;
726  if (jetCont->GetRhoParameter()) rho = jetCont->GetRhoVal();
727  if(fDebug) Printf("rho = %f",rho);
728 
729  // get MC particle container in case running embedding, to match
730  // reconstructed and MC-level trigger tracks
731  AliParticleContainer *partCont = 0x0;
733 // AliParticleContainer *partCont = 0x0;
734 // TIter next(&fParticleCollArray);
735 // while ((partCont = static_cast<AliParticleContainer*>(next()))) {
736 // TString groupname = partCont->GetName();
737 // Printf("particle name = %s",groupname.Data());
738 // }
739 
740  // Choose trigger track
741  Int_t nT=0;
742  TList ParticleList;
743  Double_t minT=0;
744  Double_t maxT=0;
745  Int_t number=0;
746  Double_t dice=fRandom->Uniform(0,1);
747  Bool_t isSignal = kFALSE;
748  if(dice>fFrac){
749  minT=fTTLowRef;
750  maxT=fTTUpRef;
751  }
752  if(dice<=fFrac){
753  isSignal = kTRUE;
754  minT=fTTLowSig;
755  maxT=fTTUpSig;
756  }
757  nT=SelectTrigger(&ParticleList,minT,maxT,number);
758  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);
759  if(nT<0) return;
760 
761  if(dice>fFrac) fh1TrigRef->Fill(number);
762  if(dice<=fFrac)fh1TrigSig->Fill(number);
763 
764 
765  // particle loop -
766  for(Int_t tt=0;tt<ParticleList.GetEntries();tt++){
767  // histogram
768  //if(fHardest==0||fHardest==1){if(tt!=nT) continue;}
769  if(tt!=nT) continue;
770  AliVParticle *partback = (AliVParticle*)ParticleList.At(tt);
771  if(!partback) continue;
772  if(fDebug) Printf("trigger particle pt = %f \teta = %f \t phi = %f",partback->Pt(),partback->Eta(),partback->Phi());
773  // if(partback->Pt()<8) continue;
774 
775  Int_t injet4=0;
776  Int_t injet=0;
777 
778  fh2Ntriggers->Fill(fCent,partback->Pt());
779  Double_t phiBinT = RelativePhi(partback->Phi(),fEPV0);
780  if(fCent<20.) fh2RPTC20->Fill(TMath::Abs(phiBinT),partback->Pt());
781  if(fCent<10.) fh2RPTC10->Fill(TMath::Abs(phiBinT),partback->Pt());
782 
783  Double_t etabig=0;
784  Double_t ptbig=0;
785  Double_t areabig=0;
786  Double_t phibig=0.;
787  // Double_t areasmall=0;
788 
789  TString histname;
790  TString groupname;
791  groupname = jetCont->GetName();
792  UInt_t count = 0;
793  for(auto jetbig : jetCont->accepted()) {
794  if (!jetbig) continue;
795  count++;
796  ptbig = jetbig->Pt();
797  etabig = jetbig->Eta();
798  phibig = jetbig->Phi();
799  if(ptbig==0) continue;
800  Double_t phiBin = RelativePhi(phibig,fEPV0); //relative phi between jet and ev. plane
801  areabig = jetbig->Area();
802  Double_t ptcorr=ptbig-rho*areabig;
803  //JJJ - perhaps should change eta selection if implemented in jet container
804  if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
805  if(areabig>=0.07) injet=injet+1;
806  if(areabig>=0.4) injet4=injet4+1;
807  Double_t dphi=RelativePhi(partback->Phi(),phibig);
808  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);
809 
810  // do azimuthal correlation analysis
811  // dPhi between -0.5 < dPhi < 1.5
812  Double_t dPhiShift=phibig-partback->Phi();
813  if(dPhiShift>2*TMath::Pi()) dPhiShift -= 2*TMath::Pi();
814  if(dPhiShift<-2*TMath::Pi()) dPhiShift += 2*TMath::Pi();
815  if(dPhiShift<-0.5*TMath::Pi()) dPhiShift += 2*TMath::Pi();
816  if(dPhiShift>1.5*TMath::Pi()) dPhiShift -= 2*TMath::Pi();
817  if(isSignal) fhDphiPtSig->Fill(dPhiShift,ptcorr);
818  else fhDphiPtRef->Fill(dPhiShift,ptcorr);
819 
820  // selection on relative phi
821  if(fJetHadronDeltaPhi>0. &&
822  TMath::Abs(dphi)<TMath::Pi()-fJetHadronDeltaPhi) continue;
823 
824  if(fCent<10.) fh2RPJetsC10->Fill(TMath::Abs(phiBin), ptcorr);
825  if(fCent<20.) fh2RPJetsC20->Fill(TMath::Abs(phiBin), ptcorr);
826 
827  Float_t phitt=partback->Phi();
828  if(phitt<0)phitt+=TMath::Pi()*2.;
829  Int_t phiBintt = GetPhiBin(phitt-fEPV0);
830 
831  Double_t fillspec[] = {fCent,areabig,ptcorr,partback->Pt(),dPhiShift, static_cast<Double_t>(phiBintt)};
832  fHJetSpec->Fill(fillspec);
833 
835  //
836  // embedding for recoil jets
837  // get MC info
838  //
839  Double_t ptTTMC = 0;
840  Double_t phiTTMC = 0;
841  Int_t TTmatched = 0;
842  for(auto partMC : partCont->accepted()) {
843  Int_t labtr = partback->GetLabel();
844  Int_t labpa = partMC->GetLabel();
845  if(labtr==labpa) {
846  ptTTMC = partMC->Pt();
847  phiTTMC = partMC->Phi();
848  TTmatched++;
849  break;
850  }
851  }
852  if(TTmatched!=1) continue;
853  Double_t ptTTreco = partback->Pt();
854  Double_t phiTTreco = partback->Phi();
855  if(fDebug) Printf("found corresponding truth-level particle, pt reco = %f pt MC = %f, phi reco = %f phi MC = %f",ptTTreco,ptTTMC,phiTTreco,phiTTMC);
856  fhTTPtDetMatchedToPart->Fill(ptTTreco,ptTTMC);
857  fhTTPhiDetMatchedToPart->Fill(phiTTreco,phiTTMC);
858 
859  auto jet2 = jetbig->ClosestJet();
860  if(!jet2) {
861  //Printf("jet 2 cant be found");
862  continue;}
863  Double_t ptJet2 = jet2->Pt();
864  fhPtDetRecoil->Fill(ptJet2);
865  auto jet3 = jet2->ClosestJet();
866  if(!jet3) {
867  //Printf("jet3 can't be found");
868  continue;
869  }
870  fhPtDetMatchedToPartRecoil->Fill(ptJet2);
871  Double_t ptJet3 = jet3->Pt();
872  Double_t phiJet3 = jet3->Phi();
873 
874 
875  Double_t dPhiPart=phiJet3-phiTTMC;
876  if(dPhiPart>2*TMath::Pi()) dPhiPart -= 2*TMath::Pi();
877  if(dPhiPart<-2*TMath::Pi()) dPhiPart += 2*TMath::Pi();
878  if(dPhiPart<-0.5*TMath::Pi()) dPhiPart += 2*TMath::Pi();
879  if(dPhiPart>1.5*TMath::Pi()) dPhiPart -= 2*TMath::Pi();
880 
881  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);
882 
883  fhPtDetPartRecoil->Fill(ptJet2,ptJet3);
884  Double_t fraction = jetCont->GetFractionSharedPt(jetbig);
885  if(fraction < fMinFractionSharedPt) continue;
886 
887  Double_t residual = (ptcorr - ptJet3) / ptJet3;
888  Double_t residualDphi = (dPhiShift - dPhiPart) / dPhiPart;
889 
890  fhPtHybrDetRecoil->Fill(ptbig,ptJet2);
891  fhPtHybrPartRecoil->Fill(ptbig,ptJet3);
892  fhPtHybrPartCorRecoil->Fill(ptcorr,ptJet3);
893  fhDPhiHybrPartCorRecoil->Fill(dPhiShift,dPhiPart);
894  fhResidualRecoil->Fill(residual);
895  fhPtResidualRecoil->Fill(ptJet3,residual);
896  fhDphiResidualRecoil->Fill(residualDphi);
897  fhDphiphiResidualRecoil->Fill(dPhiPart,residualDphi);
898 
899  if(fFillRecoilTree) {
900  fTreeVarsRecoil[0] = fCent;
901  fTreeVarsRecoil[1] = partback->Pt();
902  fTreeVarsRecoil[2] = ptbig;
903  fTreeVarsRecoil[3] = areabig;
904  fTreeVarsRecoil[4] = ptcorr;
905  fTreeVarsRecoil[5] = phibig;
906  fTreeVarsRecoil[6] = dPhiShift;
907  fTreeVarsRecoil[7] = ptJet3;
908  fTreeVarsRecoil[8] = phiJet3;
909  fTreeVarsRecoil[9] = dPhiPart;
911  fTreeEmbRecoil->Fill();
912  }
913  }
914  }
915  }
916 }
917 
919 
923  if(fDebug) Printf("particle container 0 entries = %i \t1 entries = %i\t 2 entries = %i",partCont0->GetNParticles(),partCont1->GetNParticles(),partCont2->GetNParticles());
927 
928  if(!jetCont || !jetContPart || !jetContTrue)
929  {
930  AliError(Form("jet container not found - check name %s",fJetContName.Data()));
931  TIter next(&fJetCollArray);
932  while ((jetCont = static_cast<AliJetContainer*>(next())))
933  AliError(Form("%s",jetCont->GetName()));
934  AliFatal("Exit...");
935  return;
936  }
937 
938  if(fDebug) {
939  Printf("n particle jets = %i",jetContPart->GetNJets());
940  Printf("n reco jets = %i",jetCont->GetNJets());
941  Printf("n PYTHIA jets = %i",jetContTrue->GetNJets());
942  }
943 
944  // Background
945  Double_t rho = 0;
946  if (jetCont->GetRhoParameter()) rho = jetCont->GetRhoVal();
947 
948  // PYTHIA event weight
949  // note - not used
950  // AliGenPythiaEventHeader *pyHeader = 0x0; //!<! Pythia header of the current external event
951  // AliAODEvent *ev = dynamic_cast<AliAODEvent*>(InputEvent());
952  // AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(ev->FindListObject(AliAODMCHeader::StdBranchName()));
953  // if (aodMCH) {
954  // for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
955  // pyHeader= dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
956  // if (pyHeader) break;
957  // }
958  // }
959  //
960  // Double_t pythiaCrossSection = 0;
961  // Double_t pythiaTrials = 0;
962  // Double_t pythiaWeight = 0;
963  // if (pyHeader)
964  // {
965  // if(fDebug) Printf("have pythia header - get weight");
966  // pythiaCrossSection = pyHeader->GetXsection();
967  // pythiaTrials = pyHeader->Trials();
968  // //fPythiaPtHard = fPythiaHeader->GetPtHard();
969  // pythiaWeight = pythiaCrossSection / pythiaTrials;
970  // }
971  // if(fDebug) Printf("pythia weight = %f",pythiaWeight);
972 
973  for(auto jet1 : jetCont->accepted()) { // loop over hybrid jets
974 
975  Double_t ptJet1 = jet1->Pt();
976  Double_t phiJet1 = jet1->Phi();
977  auto jet2 = jet1->ClosestJet();
978  if(!jet2) {
979  //Printf("jet 2 cant be found");
980  continue;}
981  Double_t ptJet2 = jet2->Pt();
982  fhPtDet->Fill(ptJet2);
983  auto jet3 = jet2->ClosestJet();
984  if(!jet3) {
985  //Printf("jet3 can't be found");
986  continue;
987  }
988  fhPtDetMatchedToPart->Fill(ptJet2);
989  Double_t ptJet3 = jet3->Pt();
990  Double_t phiJet3 = jet3->Phi();
991 
992  if(phiJet3>2*TMath::Pi()) phiJet3 -= 2*TMath::Pi();
993  if(phiJet3<-2*TMath::Pi()) phiJet3 += 2*TMath::Pi();
994  if(phiJet3<-0.5*TMath::Pi()) phiJet3 += 2*TMath::Pi();
995  if(phiJet3>1.5*TMath::Pi()) phiJet3 -= 2*TMath::Pi();
996 
997  if(phiJet1>2*TMath::Pi()) phiJet1 -= 2*TMath::Pi();
998  if(phiJet1<-2*TMath::Pi()) phiJet1 += 2*TMath::Pi();
999  if(phiJet1<-0.5*TMath::Pi()) phiJet1 += 2*TMath::Pi();
1000  if(phiJet1>1.5*TMath::Pi()) phiJet1 -= 2*TMath::Pi();
1001 
1002  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);
1003 
1004  fhPtDetPart->Fill(ptJet2,ptJet3);
1005  Double_t fraction = jetCont->GetFractionSharedPt(jet1);
1006  if(fraction < fMinFractionSharedPt) continue;
1007 
1008  Double_t area = jet1->Area();
1009  Double_t ptCorr = ptJet1-rho*area;
1010 
1011  Double_t residual = (ptCorr - ptJet3) / ptJet3;
1012  Double_t residualPhi = (phiJet1 - phiJet3) / phiJet3;
1013 
1014  fhPtHybrDet->Fill(ptJet1,ptJet2);
1015  fhPtHybrPart->Fill(ptJet1,ptJet3);
1016  fhPtHybrPartCor->Fill(ptCorr,ptJet3);
1017  fhPhiHybrPartCor->Fill(phiJet1,phiJet3);
1018 
1019  fhResidual->Fill(residual);
1020  fhPtResidual->Fill(ptJet3,residual);
1021  fhPhiResidual->Fill(residualPhi);
1022  fhPhiPhiResidual->Fill(phiJet3,residualPhi);
1023 
1024  if(fFillInclusiveTree) {
1026  fTreeVarsInclusive[1] = ptJet1;
1027  fTreeVarsInclusive[2] = area;
1028  fTreeVarsInclusive[3] = ptCorr;
1029  fTreeVarsInclusive[4] = phiJet1;
1030  fTreeVarsInclusive[5] = ptJet3;
1031  fTreeVarsInclusive[6] = phiJet3;
1033  fTreeEmbInclusive->Fill();
1034  }
1035  }
1036 // for(auto jettrue : jetContPart->accepted()) {
1037 // // jettrue
1038 // Double_t ptTrue = jettrue->Pt();
1039 // Double_t phiTrue = jettrue->Phi();
1040 // Double_t etaTrue= jettrue->Eta();
1041 // auto jetmatched = jettrue->ClosestJet();
1042 // if(!jetmatched) continue;
1043 // Double_t ptMatched = jetmatched->Pt();
1044 // Double_t phiMatched = jetmatched->Phi();
1045 // Double_t etaMatched = jetmatched->Eta();
1046 // Double_t residual = (ptMatched - ptTrue) / ptTrue;
1047 // fhPtHybrTrue->Fill(ptTrue,ptMatched);
1048 // fhResidual->Fill(residual);
1049 // fhPtResidual->Fill(ptTrue,residual);
1050 // }
1051 }
1052 
1053 
1054 
1060 {
1061  TString histname;
1062  TString groupname;
1063  AliJetContainer* jetCont = 0;
1064  TIter next(&fJetCollArray);
1065  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
1066  groupname = jetCont->GetName();
1067  UInt_t count = 0;
1068  for(auto jet : jetCont->accepted()) {
1069  if (!jet) continue;
1070  count++;
1071 
1072  histname = TString::Format("%s/histJetPt_%d", groupname.Data(), fCentBin);
1073  fHistManager.FillTH1(histname, jet->Pt());
1074 
1075  histname = TString::Format("%s/histJetArea_%d", groupname.Data(), fCentBin);
1076  fHistManager.FillTH1(histname, jet->Area());
1077 
1078  histname = TString::Format("%s/histJetPhi_%d", groupname.Data(), fCentBin);
1079  fHistManager.FillTH1(histname, jet->Phi());
1080 
1081  histname = TString::Format("%s/histJetEta_%d", groupname.Data(), fCentBin);
1082  fHistManager.FillTH1(histname, jet->Eta());
1083 
1084  if (jetCont->GetRhoParameter()) {
1085  histname = TString::Format("%s/histJetCorrPt_%d", groupname.Data(), fCentBin);
1086  fHistManager.FillTH1(histname, jet->Pt() - jetCont->GetRhoVal() * jet->Area());
1087  }
1088  }
1089  histname = TString::Format("%s/histNJets_%d", groupname.Data(), fCentBin);
1090  fHistManager.FillTH1(histname, count);
1091  }
1092 }
1093 
1099 {
1101 
1102  TString histname;
1103  TString groupname;
1104  UInt_t sumAcceptedTracks = 0;
1105  AliParticleContainer* partCont = 0;
1106  TIter next(&fParticleCollArray);
1107  while ((partCont = static_cast<AliParticleContainer*>(next()))) {
1108  groupname = partCont->GetName();
1109  UInt_t count = 0;
1110  for(auto part : partCont->accepted()) {
1111  if (!part) continue;
1112  count++;
1113 
1114  histname = TString::Format("%s/histTrackPt_%d", groupname.Data(), fCentBin);
1115  fHistManager.FillTH1(histname, part->Pt());
1116 
1117  histname = TString::Format("%s/histTrackPhi_%d", groupname.Data(), fCentBin);
1118  fHistManager.FillTH1(histname, part->Phi());
1119 
1120  histname = TString::Format("%s/histTrackEta_%d", groupname.Data(), fCentBin);
1121  fHistManager.FillTH1(histname, part->Eta());
1122 
1123  if (partCont->GetLoadedClass()->InheritsFrom("AliVTrack")) {
1124  const AliVTrack* track = static_cast<const AliVTrack*>(part);
1125 
1126  histname = TString::Format("%s/fHistDeltaEtaPt_%d", groupname.Data(), fCentBin);
1127  fHistManager.FillTH1(histname, track->Pt(), track->Eta() - track->GetTrackEtaOnEMCal());
1128 
1129  histname = TString::Format("%s/fHistDeltaPhiPt_%d", groupname.Data(), fCentBin);
1130  fHistManager.FillTH1(histname, track->Pt(), track->Phi() - track->GetTrackPhiOnEMCal());
1131 
1132  histname = TString::Format("%s/fHistDeltaPtvsPt_%d", groupname.Data(), fCentBin);
1133  fHistManager.FillTH1(histname, track->Pt(), track->Pt() - track->GetTrackPtOnEMCal());
1134 
1135  if (clusCont) {
1136  Int_t iCluster = track->GetEMCALcluster();
1137  if (iCluster >= 0) {
1138  AliVCluster* cluster = clusCont->GetAcceptCluster(iCluster);
1139  if (cluster) {
1140  histname = TString::Format("%s/fHistEoverPvsP_%d", groupname.Data(), fCentBin);
1141  fHistManager.FillTH2(histname, track->P(), cluster->GetNonLinCorrEnergy() / track->P());
1142  }
1143  }
1144  }
1145  }
1146  }
1147  sumAcceptedTracks += count;
1148 
1149  histname = TString::Format("%s/histNTracks_%d", groupname.Data(), fCentBin);
1150  fHistManager.FillTH1(histname, count);
1151  }
1152 
1153  histname = "fHistSumNTracks";
1154  fHistManager.FillTH1(histname, sumAcceptedTracks);
1155 }
1156 
1162 {
1163  TString histname;
1164  TString groupname;
1165  UInt_t sumAcceptedClusters = 0;
1166  AliClusterContainer* clusCont = 0;
1167  TIter next(&fClusterCollArray);
1168  while ((clusCont = static_cast<AliClusterContainer*>(next()))) {
1169  groupname = clusCont->GetName();
1170 
1171  for(auto cluster : clusCont->all()) {
1172  if (!cluster) continue;
1173 
1174  if (cluster->GetIsExotic()) {
1175  histname = TString::Format("%s/histClusterEnergyExotic_%d", groupname.Data(), fCentBin);
1176  fHistManager.FillTH1(histname, cluster->E());
1177  }
1178  }
1179 
1180  UInt_t count = 0;
1181  for(auto cluster : clusCont->accepted()) {
1182  if (!cluster) continue;
1183  count++;
1184 
1185  AliTLorentzVector nPart;
1186  cluster->GetMomentum(nPart, fVertex);
1187 
1188  histname = TString::Format("%s/histClusterEnergy_%d", groupname.Data(), fCentBin);
1189  fHistManager.FillTH1(histname, cluster->E());
1190 
1191  histname = TString::Format("%s/histClusterNonLinCorrEnergy_%d", groupname.Data(), fCentBin);
1192  fHistManager.FillTH1(histname, cluster->GetNonLinCorrEnergy());
1193 
1194  histname = TString::Format("%s/histClusterHadCorrEnergy_%d", groupname.Data(), fCentBin);
1195  fHistManager.FillTH1(histname, cluster->GetHadCorrEnergy());
1196 
1197  histname = TString::Format("%s/histClusterPhi_%d", groupname.Data(), fCentBin);
1198  fHistManager.FillTH1(histname, nPart.Phi_0_2pi());
1199 
1200  histname = TString::Format("%s/histClusterEta_%d", groupname.Data(), fCentBin);
1201  fHistManager.FillTH1(histname, nPart.Eta());
1202  }
1203  sumAcceptedClusters += count;
1204 
1205  histname = TString::Format("%s/histNClusters_%d", groupname.Data(), fCentBin);
1206  fHistManager.FillTH1(histname, count);
1207  }
1208 
1209  histname = "fHistSumNClusters";
1210  fHistManager.FillTH1(histname, sumAcceptedClusters);
1211 }
1212 
1218 {
1219  if (!fCaloCells) return;
1220 
1221  TString histname;
1222 
1223  const Short_t ncells = fCaloCells->GetNumberOfCells();
1224 
1225  histname = TString::Format("%s/histNCells_%d", fCaloCellsName.Data(), fCentBin);
1226  fHistManager.FillTH1(histname, ncells);
1227 
1228  histname = TString::Format("%s/histCellEnergy_%d", fCaloCellsName.Data(), fCentBin);
1229  for (Short_t pos = 0; pos < ncells; pos++) {
1230  Double_t amp = fCaloCells->GetAmplitude(pos);
1231 
1232  fHistManager.FillTH1(histname, amp);
1233  }
1234 }
1235 
1241 {
1243 }
1244 
1253 {
1254  return kTRUE;
1255 }
1256 
1261 {
1262 }
1263 
1265 
1266  Int_t index=-1;
1267  Int_t triggers[100];
1268 
1269  for(Int_t cr=0;cr<100;cr++) triggers[cr]=-1;
1270 
1271  Int_t im=0;
1272 
1273  TString groupname = "";
1274  AliParticleContainer* partCont = 0x0;
1276  else partCont = GetParticleContainer(0);
1277  groupname = partCont->GetName();
1278  UInt_t iCount = 0;
1279  for(auto part : partCont->accepted()) {
1280  if (!part) continue;
1281  list->Add(part);
1282  iCount++;
1283  if(part->Pt()>=minT && part->Pt()<maxT){
1284  triggers[im]=iCount-1;
1285  im=im+1;
1286  }
1287  }
1288  number=im;
1289  Int_t rd=0;
1290  if(im>0) rd=fRandom->Integer(im);
1291  index=triggers[rd];
1292  return index;
1293 }
1294 
1296 
1297  // get relative DeltaPhi from -pi < DeltaPhi < pi
1298 
1299  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1300  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1301  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1302  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1303  double dphi = mphi-vphi;
1304  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1305  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1306  return dphi;//dphi in [-Pi, Pi]
1307 }
1308 
1310 {
1311  Int_t phibin=-1;
1312  if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1313  Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1314  phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1315  if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1316  return phibin;
1317 }
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
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
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