AliPhysics  cc1c0ba (cc1c0ba)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskEmcalVsPhos.cxx
Go to the documentation of this file.
1 /**********************************************************************************
2 * Copyright (C) 2016, Copyright Holders of the ALICE Collaboration *
3 * All rights reserved. *
4 * *
5 * Redistribution and use in source and binary forms, with or without *
6 * modification, are permitted provided that the following conditions are met: *
7 * * Redistributions of source code must retain the above copyright *
8 * notice, this list of conditions and the following disclaimer. *
9 * * Redistributions in binary form must reproduce the above copyright *
10 * notice, this list of conditions and the following disclaimer in the *
11 * documentation and/or other materials provided with the distribution. *
12 * * Neither the name of the <organization> nor the *
13 * names of its contributors may be used to endorse or promote products *
14 * derived from this software without specific prior written permission. *
15 * *
16 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND *
17 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED *
18 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE *
19 * DISCLAIMED. IN NO EVENT SHALL ALICE COLLABORATION BE LIABLE FOR ANY *
20 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES *
21 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; *
22 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND *
23 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT *
24 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
25 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
26 * *********************************************************************************/
27 
28 #include <vector>
29 
30 #include <TClonesArray.h>
31 #include <TH1F.h>
32 #include <TH2F.h>
33 #include <TH3.h>
34 #include <TList.h>
35 #include <THnSparse.h>
36 
37 #include <AliVCluster.h>
38 #include <AliVParticle.h>
39 #include <AliLog.h>
40 
41 #include "AliTLorentzVector.h"
42 #include "AliEmcalJet.h"
43 #include "AliRhoParameter.h"
44 #include "AliJetContainer.h"
45 #include "AliParticleContainer.h"
46 #include "AliClusterContainer.h"
47 #include "AliEMCALGeometry.h"
48 #include "AliPHOSGeometry.h"
49 #include "AliOADBContainer.h"
50 #include "AliMCAnalysisUtils.h"
51 
53 
57 
63  fPlotClusterHistograms(kTRUE),
64  fPlotNeutralJets(kFALSE),
65  fPlotClustersInJets(kFALSE),
66  fPlotCellHistograms(kTRUE),
67  fPlotClusWithoutNonLinCorr(kFALSE),
68  fPlotExotics(kFALSE),
69  fPlotStandardClusterTHnSparse(kTRUE),
70  fPlotNearestNeighborDistribution(kFALSE),
71  fPlotClusterCone(kFALSE),
72  fPlotCaloCentrality(kFALSE),
73  fPlotFineGrainedEtaPhi(kFALSE),
74  fPlotEvenOddEta(kFALSE),
75  fPlotCellSMDensity(kFALSE),
76  fExcludeRejectedCells(kFALSE),
77  fPlotFineGrainedCentrality(kFALSE),
78  fPlotJetPerformance(kFALSE),
79  fMaxPt(200),
80  fNCentHistBins(0),
81  fCentHistBins(0),
82  fNPtHistBins(0),
83  fPtHistBins(0),
84  fNM02HistBins(0),
85  fM02HistBins(0),
86  fUseAliEventCuts(kTRUE),
87  fEventCuts(0),
88  fEventCutList(0),
89  fUseManualEventCuts(kFALSE),
90  fGeneratorLevel(0),
91  fPHOSGeo(nullptr),
92  fHistManager()
93 {
95 }
96 
103  AliAnalysisTaskEmcalJet(name, kTRUE),
104  fPlotClusterHistograms(kTRUE),
105  fPlotNeutralJets(kFALSE),
106  fPlotClustersInJets(kFALSE),
107  fPlotCellHistograms(kTRUE),
108  fPlotClusWithoutNonLinCorr(kFALSE),
109  fPlotExotics(kFALSE),
110  fPlotStandardClusterTHnSparse(kTRUE),
111  fPlotNearestNeighborDistribution(kFALSE),
112  fPlotClusterCone(kFALSE),
113  fPlotCaloCentrality(kFALSE),
114  fPlotFineGrainedEtaPhi(kFALSE),
115  fPlotEvenOddEta(kFALSE),
116  fPlotCellSMDensity(kFALSE),
117  fExcludeRejectedCells(kFALSE),
118  fPlotFineGrainedCentrality(kFALSE),
119  fPlotJetPerformance(kFALSE),
120  fMaxPt(200),
121  fNCentHistBins(0),
122  fCentHistBins(0),
123  fNPtHistBins(0),
124  fPtHistBins(0),
125  fNM02HistBins(0),
126  fM02HistBins(0),
127  fUseAliEventCuts(kTRUE),
128  fEventCuts(0),
129  fEventCutList(0),
130  fUseManualEventCuts(kFALSE),
131  fGeneratorLevel(0),
132  fPHOSGeo(nullptr),
133  fHistManager(name)
134 {
136 }
137 
142 {
143 }
144 
149 {
150  fNCentHistBins = 4;
152  fCentHistBins[0] = 0;
153  fCentHistBins[1] = 10;
154  fCentHistBins[2] = 30;
155  fCentHistBins[3] = 50;
156  fCentHistBins[4] = 90;
157 
158  fNPtHistBins = 82;
161  GenerateFixedBinArray(7, 0.3, 1, fPtHistBins+6);
162  GenerateFixedBinArray(10, 1, 3, fPtHistBins+13);
163  GenerateFixedBinArray(14, 3, 10, fPtHistBins+23);
164  GenerateFixedBinArray(10, 10, 20, fPtHistBins+37);
165  GenerateFixedBinArray(15, 20, 50, fPtHistBins+47);
166  GenerateFixedBinArray(20, 50, 150, fPtHistBins+62);
167 
168  fNM02HistBins = 81;
170  GenerateFixedBinArray(35, 0, 0.7, fM02HistBins);
171  GenerateFixedBinArray(6, 0.7, 1., fM02HistBins+35);
172  GenerateFixedBinArray(20, 1., 3., fM02HistBins+41);
173  GenerateFixedBinArray(10, 3., 5., fM02HistBins+61);
174  GenerateFixedBinArray(10, 5., 10., fM02HistBins+71);
175 }
176 
182 {
184 
185  fGeneratorLevel = GetMCParticleContainer("mcparticles");
186 
188 
189  TIter next(fHistManager.GetListOfHistograms());
190  TObject* obj = 0;
191  while ((obj = next())) {
192  fOutput->Add(obj);
193  }
194 
195  // Intialize AliEventCuts
196  if (fUseAliEventCuts) {
197  fEventCutList = new TList();
198  fEventCutList ->SetOwner();
199  fEventCutList ->SetName("EventCutOutput");
200 
201  fEventCuts.OverrideAutomaticTriggerSelection(fOffTrigger);
202  if(fUseManualEventCuts==1)
203  {
204  fEventCuts.SetManualMode();
205  // Configure manual settings here
206  // ...
207  }
208  fEventCuts.AddQAplotsToList(fEventCutList);
209  fOutput->Add(fEventCutList);
210  }
211 
212  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
213 }
214 
215 /*
216  * This function allocates the histograms for the calorimeter performance study.
217  */
219 {
220 
223  }
224 
225  if (fPlotCellHistograms) {
227  }
228 
229  if (fPlotNeutralJets) {
231  }
232 
233  if (fPlotClustersInJets) {
235  }
236 
237  if (fPlotJetPerformance) {
239  }
240 
241 }
242 
243 /*
244  * This function allocates the histograms for the calorimeter performance study.
245  */
247 {
248  TString histname;
249  TString htitle;
250 
251  Double_t* clusType = new Double_t[3+1];
252  GenerateFixedBinArray(3, -0.5, 2.5, clusType);
253  const Int_t nRejBins = 32;
254  Double_t* rejReasonBins = new Double_t[nRejBins+1];
255  GenerateFixedBinArray(nRejBins, 0, nRejBins, rejReasonBins);
256  const Int_t nExBins = 200;
257  Double_t* exBins = new Double_t[nExBins+1];
258  GenerateFixedBinArray(nExBins, 0, 1, exBins);
259 
260  AliEmcalContainer* cont = 0;
261  TIter nextClusColl(&fClusterCollArray);
262  while ((cont = static_cast<AliEmcalContainer*>(nextClusColl()))) {
263 
264  // rejection reason plot, to make efficiency correction
265  histname = TString::Format("%s/hClusterRejectionReasonEMCal", cont->GetArrayName().Data());
266  htitle = histname + ";Rejection reason;#it{E}_{clus} (GeV/)";
267  TH2* hist = fHistManager.CreateTH2(histname.Data(), htitle.Data(), nRejBins, rejReasonBins, fNPtHistBins, fPtHistBins);
268  SetRejectionReasonLabels(hist->GetXaxis());
269 
270  histname = TString::Format("%s/hClusterRejectionReasonPHOS", cont->GetArrayName().Data());
271  htitle = histname + ";Rejection reason;#it{E}_{clus} (GeV/)";
272  TH2* histPhos = fHistManager.CreateTH2(histname.Data(), htitle.Data(), nRejBins, rejReasonBins, fNPtHistBins, fPtHistBins);
273  SetRejectionReasonLabels(histPhos->GetXaxis());
274 
275  histname = TString::Format("%s/hClusterRejectionReasonMC", cont->GetArrayName().Data());
276  htitle = histname + ";Rejection reason;#it{E}_{clus} (GeV/)";
277  TH2* histMC = fHistManager.CreateTH2(histname.Data(), htitle.Data(), nRejBins, rejReasonBins, fNPtHistBins, fPtHistBins);
278  SetRejectionReasonLabels(histMC->GetXaxis());
279 
280  // plot by SM
281  const Int_t nEmcalSM = 20;
282  for (Int_t sm = 0; sm < nEmcalSM; sm++) {
283  histname = TString::Format("%s/BySM/hEmcalClusEnergy_SM%d", cont->GetArrayName().Data(), sm);
284  htitle = histname + ";#it{E}_{cluster} (GeV);counts";
285  fHistManager.CreateTH1(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins);
286  }
287 
288  for (Int_t sm = 1; sm < 5; sm++) {
289  histname = TString::Format("%s/BySM/hPhosClusEnergy_SM%d", cont->GetArrayName().Data(), sm);
290  htitle = histname + ";#it{E}_{cluster} (GeV);counts";
291  fHistManager.CreateTH1(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins);
292  }
293 
294  // Plot cluster THnSparse (centrality, eta, phi, E, E-hadcorr, has matched track, M02, Ncells, passed dispersion cut)
295  Int_t dim = 0;
296  TString title[20];
297  Int_t nbins[20] = {0};
298  Double_t min[30] = {0.};
299  Double_t max[30] = {0.};
300  Double_t *binEdges[20] = {0};
301 
303  title[dim] = "Centrality %";
305  nbins[dim] = 18;
306  min[dim] = 0;
307  max[dim] = 90;
308  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
309  }
310  else {
311  nbins[dim] = fNCentHistBins;
312  binEdges[dim] = fCentHistBins;
313  min[dim] = fCentHistBins[0];
314  max[dim] = fCentHistBins[fNCentHistBins];
315  }
316  dim++;
317  }
318 
319  title[dim] = "#eta";
320  nbins[dim] = 28;
322  nbins[dim] = 100;
323  }
324  min[dim] = -0.7;
325  max[dim] = 0.7;
326  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
327  dim++;
328 
329  title[dim] = "#phi";
330  nbins[dim] = 100;
332  nbins[dim] = 357;
333  }
334  min[dim] = 1.;
335  max[dim] = 6.;
336  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
337  dim++;
338 
339  title[dim] = "#it{E}_{clus} (GeV)";
340  nbins[dim] = fNPtHistBins;
341  binEdges[dim] = fPtHistBins;
342  min[dim] = fPtHistBins[0];
343  max[dim] = fPtHistBins[fNPtHistBins];
344  dim++;
345 
347  title[dim] = "#it{E}_{clus, hadcorr} or #it{E}_{core} (GeV)";
348  nbins[dim] = fNPtHistBins;
349  binEdges[dim] = fPtHistBins;
350  min[dim] = fPtHistBins[0];
351  max[dim] = fPtHistBins[fNPtHistBins];
352  dim++;
353 
354  title[dim] = "Matched track";
355  nbins[dim] = 2;
356  min[dim] = -0.5;
357  max[dim] = 1.5;
358  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
359  dim++;
360 
361  title[dim] = "M02";
362  nbins[dim] = fNM02HistBins;
363  binEdges[dim] = fM02HistBins;
364  min[dim] = fM02HistBins[0];
365  max[dim] = fM02HistBins[fNM02HistBins];
366  dim++;
367 
368  title[dim] = "Ncells";
369  nbins[dim] = 30;
370  min[dim] = -0.5;
371  max[dim] = 29.5;
372  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
373  dim++;
374 
375  title[dim] = "Dispersion cut";
376  nbins[dim] = 2;
377  min[dim] = -0.5;
378  max[dim] = 1.5;
379  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
380  dim++;
381 
382  if (fGeneratorLevel) {
383  title[dim] = "Particle type1";
384  nbins[dim] = 8;
385  min[dim] = -0.5;
386  max[dim] = 7.5;
387  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
388  dim++;
389 
390  title[dim] = "Particle type2";
391  nbins[dim] = 8;
392  min[dim] = -0.5;
393  max[dim] = 7.5;
394  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
395  dim++;
396  }
397  }
398 
399  if (fPlotEvenOddEta) {
400  title[dim] = "Even/odd eta";
401  nbins[dim] = 2;
402  min[dim] = -0.5;
403  max[dim] = 1.5;
404  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
405  dim++;
406  }
407 
409  title[dim] = "#DeltaR_{NN}";
410  nbins[dim] = 100;
411  min[dim] = 0.;
412  max[dim] = 1.;
413  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
414  dim++;
415  }
416 
417  if (fPlotClusterCone) {
418 
419  title[dim] = "Cone type";
420  nbins[dim] = 2;
421  min[dim] = -0.5;
422  max[dim] = 1.5;
423  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
424  dim++;
425 
426  title[dim] = "R";
427  nbins[dim] = 2;
428  min[dim] = 0;
429  max[dim] = 0.15;
430  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
431  dim++;
432 
433  title[dim] = "#it{E}_{cone} (GeV)";
434  nbins[dim] = fNPtHistBins;
435  binEdges[dim] = fPtHistBins;
436  min[dim] = fPtHistBins[0];
437  max[dim] = fPtHistBins[fNPtHistBins];
438  dim++;
439 
440  }
441 
442  if (fPlotCaloCentrality) {
443 
444  title[dim] = "#it{#rho}_{cell cone} (GeV)";
445  nbins[dim] = 100;
446  min[dim] = 0.;
447  max[dim] = 1000.;
448  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
449  dim++;
450 
451  title[dim] = "Ncells cone";
452  nbins[dim] = 100;
453  min[dim] = -0.5;
454  max[dim] = 99.5;
455  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
456  dim++;
457 
458  if (fPlotCellSMDensity) {
459  title[dim] = "#it{E}_{cell SM} (GeV)";
460  nbins[dim] = fNPtHistBins;
461  binEdges[dim] = fPtHistBins;
462  min[dim] = fPtHistBins[0];
463  max[dim] = fPtHistBins[fNPtHistBins];
464  dim++;
465 
466  title[dim] = "Ncells SM";
467  nbins[dim] = 100;
468  min[dim] = -0.5;
469  max[dim] = 999.5;
470  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
471  dim++;
472  }
473 
474  }
475 
476  TString thnname = TString::Format("%s/clusterObservables", cont->GetArrayName().Data());
477  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
478  for (Int_t i = 0; i < dim; i++) {
479  hn->GetAxis(i)->SetTitle(title[i]);
480  hn->SetBinEdges(i, binEdges[i]);
481  }
482 
483  // Plot Fcross distribution
484  if (fPlotExotics) {
485  histname = TString::Format("%s/hFcrossEMCal", cont->GetArrayName().Data());
486  htitle = histname + ";Centrality (%);Fcross;#it{E}_{clus} (GeV/)";
487  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, nExBins, exBins, fNPtHistBins, fPtHistBins);
488  }
489 
490  }
491 
492 }
493 
494 /*
495  * This function allocates the histograms for calo cells.
496  */
498 {
499  TString histname;
500  TString htitle;
501 
502  Double_t* clusType = new Double_t[3+1];
503  GenerateFixedBinArray(3, -0.5, 2.5, clusType);
504 
505  // centrality vs. cell energy vs. cell type (for all cells)
506  histname = TString::Format("Cells/hCellEnergyAll");
507  htitle = histname + ";#it{E}_{cell} (GeV);Centrality (%); Cluster type";
508  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins, 3, clusType);
509 
510  // centrality vs. cell energy vs. cell type (for cells in accepted clusters)
511  histname = TString::Format("Cells/hCellEnergyAccepted");
512  htitle = histname + ";#it{E}_{cell} (GeV);Centrality (%); Cluster type";
513  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins, 3, clusType);
514 
515  // centrality vs. cell energy vs. cell type (for leading cells in accepted clusters)
516  histname = TString::Format("Cells/hCellEnergyLeading");
517  htitle = histname + ";#it{E}_{cell} (GeV);Centrality (%); Cluster type";
518  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins, 3, clusType);
519 
520  // plot cell patches by SM
521  const Int_t nEmcalSM = 20;
522  for (Int_t sm = 0; sm < nEmcalSM; sm++) {
523  histname = TString::Format("Cells/BySM/hEmcalPatchEnergy_SM%d", sm);
524  htitle = histname + ";#it{E}_{cell patch} (GeV);Centrality (%)";
525  fHistManager.CreateTH2(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins);
526  }
527 
528  for (Int_t sm = 1; sm < 5; sm++) {
529  histname = TString::Format("Cells/BySM/hPhosPatchEnergy_SM%d", sm);
530  htitle = histname + ";#it{E}_{cell patch} (GeV);Centrality (%)";
531  fHistManager.CreateTH2(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins);
532  }
533 
534 }
535 
536 /*
537  * This function allocates the histograms for neutral jets.
538  */
540 {
541  AliJetContainer* jets = 0;
542  TIter nextJetColl(&fJetCollArray);
543  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
544 
545  TString axisTitle[30]= {""};
546  Int_t nbinsJet[30] = {0};
547  Double_t minJet[30] = {0.};
548  Double_t maxJet[30] = {0.};
549  Double_t *binEdgesJet[20] = {0};
550  Int_t dimJet = 0;
551 
552  if (fForceBeamType != kpp) {
553  axisTitle[dimJet] = "Centrality (%)";
554  nbinsJet[dimJet] = fNCentHistBins;
555  binEdgesJet[dimJet] = fCentHistBins;
556  minJet[dimJet] = fCentHistBins[0];
557  maxJet[dimJet] = fCentHistBins[fNCentHistBins];
558  dimJet++;
559  }
560 
561  axisTitle[dimJet] = "#eta_{jet}";
562  nbinsJet[dimJet] = 28;
563  minJet[dimJet] = -0.7;
564  maxJet[dimJet] = 0.7;
565  binEdgesJet[dimJet] = GenerateFixedBinArray(nbinsJet[dimJet], minJet[dimJet], maxJet[dimJet]);
566  dimJet++;
567 
568  axisTitle[dimJet] = "#phi_{jet} (rad)";
569  nbinsJet[dimJet] = 100;
570  minJet[dimJet] = 1.;
571  maxJet[dimJet] = 6.;
572  binEdgesJet[dimJet] = GenerateFixedBinArray(nbinsJet[dimJet], minJet[dimJet], maxJet[dimJet]);
573  dimJet++;
574 
575  axisTitle[dimJet] = "#it{E}_{T} (GeV)";
576  nbinsJet[dimJet] = fNPtHistBins;
577  binEdgesJet[dimJet] = fPtHistBins;
578  minJet[dimJet] = fPtHistBins[0];
579  maxJet[dimJet] = fPtHistBins[fNPtHistBins];
580  dimJet++;
581 
582  axisTitle[dimJet] = "#rho (GeV/#it{c})";
583  nbinsJet[dimJet] = 100;
584  minJet[dimJet] = 0.;
585  maxJet[dimJet] = 1000.;
586  binEdgesJet[dimJet] = GenerateFixedBinArray(nbinsJet[dimJet], minJet[dimJet], maxJet[dimJet]);
587  dimJet++;
588 
589  axisTitle[dimJet] = "N_{clusters}";
590  nbinsJet[dimJet] = 20;
591  minJet[dimJet] = -0.5;
592  maxJet[dimJet] = 19.5;
593  binEdgesJet[dimJet] = GenerateFixedBinArray(nbinsJet[dimJet], minJet[dimJet], maxJet[dimJet]);
594  dimJet++;
595 
596  axisTitle[dimJet] = "#it{E}_{T}, acc clus within R (GeV)";
597  nbinsJet[dimJet] = fNPtHistBins;
598  binEdgesJet[dimJet] = fPtHistBins;
599  minJet[dimJet] = fPtHistBins[0];
600  maxJet[dimJet] = fPtHistBins[fNPtHistBins];
601  dimJet++;
602 
603  axisTitle[dimJet] = "#it{E}_{T}, acc cell within R (GeV)";
604  nbinsJet[dimJet] = fNPtHistBins;
605  binEdgesJet[dimJet] = fPtHistBins;
606  minJet[dimJet] = fPtHistBins[0];
607  maxJet[dimJet] = fPtHistBins[fNPtHistBins];
608  dimJet++;
609 
610  TString thnname = TString::Format("%s/hClusterJetObservables", jets->GetArrayName().Data());
611  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dimJet, nbinsJet, minJet, maxJet);
612  for (Int_t i = 0; i < dimJet; i++) {
613  hn->GetAxis(i)->SetTitle(axisTitle[i]);
614  hn->SetBinEdges(i, binEdgesJet[i]);
615  }
616 
617  }
618 }
619 
620 /*
621  * This function allocates the histograms for clusters within jets.
622  */
624 {
625  TString histname;
626  TString htitle;
627  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
628 
629  AliJetContainer* jets = 0;
630  TIter nextJetColl(&fJetCollArray);
631  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
632 
633  // Plot cluster spectra within jets
634  // (centrality, cluster energy, jet pT, jet eta, jet phi)
635  Int_t dim = 0;
636  TString title[20];
637  Int_t nbins[20] = {0};
638  Double_t min[30] = {0.};
639  Double_t max[30] = {0.};
640  Double_t *binEdges[20] = {0};
641 
643  title[dim] = "Centrality %";
644  nbins[dim] = fNCentHistBins;
645  binEdges[dim] = fCentHistBins;
646  min[dim] = fCentHistBins[0];
647  max[dim] = fCentHistBins[fNCentHistBins];
648  dim++;
649  }
650 
651  title[dim] = "#it{E}_{clus} (GeV)";
652  nbins[dim] = fNPtHistBins;
653  binEdges[dim] = fPtHistBins;
654  min[dim] = fPtHistBins[0];
655  max[dim] = fPtHistBins[fNPtHistBins];
656  dim++;
657 
658  title[dim] = "#it{p}_{T,jet}^{corr}";
659  nbins[dim] = nPtBins;
660  min[dim] = 0;
661  max[dim] = fMaxPt;
662  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
663  dim++;
664 
665  title[dim] = "#eta_{jet}";
666  nbins[dim] = 40;
667  min[dim] = -0.5;
668  max[dim] = 0.5;
669  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
670  dim++;
671 
672  title[dim] = "#phi_{jet}";
673  nbins[dim] = 200;
674  min[dim] = 1.;
675  max[dim] = 6.;
676  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
677  dim++;
678 
679  TString thnname = TString::Format("%s/hClustersInJets", jets->GetArrayName().Data());
680  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
681  for (Int_t i = 0; i < dim; i++) {
682  hn->GetAxis(i)->SetTitle(title[i]);
683  hn->SetBinEdges(i, binEdges[i]);
684  }
685 
686  // (jet type, jet pT, cluster shift)
687  histname = TString::Format("%s/hCaloJESshift", jets->GetArrayName().Data());
688  htitle = histname + ";type;#it{p}_{T}^{corr} (GeV/#it{c});#Delta_{JES}";
689  fHistManager.CreateTH3(histname.Data(), htitle.Data(), 3, -0.5, 2.5, nPtBins, 0, fMaxPt, 100, 0, 20);
690 
691  }
692 }
693 
694 /*
695  * This function allocates the histograms for the jet performance study.
696  */
698 {
699  TString histname;
700  TString htitle;
701  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
702 
703  const Int_t nParticleTypes = 8;
704  Double_t *particleTypeBins = GenerateFixedBinArray(nParticleTypes, -0.5, 7.5);
705 
706  // M02 vs. Energy vs. Particle type
707  histname = "JetPerformance/hM02VsParticleType";
708  htitle = histname + ";M02;#it{E}_{clus} (GeV); Particle type";
709  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNM02HistBins, fM02HistBins, fNPtHistBins, fPtHistBins, nParticleTypes, particleTypeBins);
710 
711  // M02 vs. Energy vs. Particle type vs. Jet pT, for particles inside jets
712  Int_t dim = 0;
713  TString title[20];
714  Int_t nbins[20] = {0};
715  Double_t min[30] = {0.};
716  Double_t max[30] = {0.};
717  Double_t *binEdges[20] = {0};
718 
719  title[dim] = "M02";
720  nbins[dim] = fNM02HistBins;
721  binEdges[dim] = fM02HistBins;
722  min[dim] = fM02HistBins[0];
723  max[dim] = fM02HistBins[fNM02HistBins];
724  dim++;
725 
726  title[dim] = "#it{E}_{clus} (GeV)";
727  nbins[dim] = fNPtHistBins;
728  binEdges[dim] = fPtHistBins;
729  min[dim] = fPtHistBins[0];
730  max[dim] = fPtHistBins[fNPtHistBins];
731  dim++;
732 
733  title[dim] = "Particle type";
734  nbins[dim] = nParticleTypes;
735  min[dim] = -0.5;
736  max[dim] = 7.5;
737  binEdges[dim] = particleTypeBins;
738  dim++;
739 
740  title[dim] = "#it{p}_{T,jet}^{corr}";
741  nbins[dim] = nPtBins;
742  min[dim] = 0;
743  max[dim] = fMaxPt;
744  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
745  dim++;
746 
747  TString thnname = "JetPerformance/hM02VsParticleTypeJets";
748  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
749  for (Int_t i = 0; i < dim; i++) {
750  hn->GetAxis(i)->SetTitle(title[i]);
751  hn->SetBinEdges(i, binEdges[i]);
752  }
753 
754  // Particle composition inside each jet -- jet pT vs. particle type vs. particle number vs. particle pT sum
755  // (One entry per jet for each particle type)
756  dim = 0;
757 
758  title[dim] = "#it{p}_{T,jet}^{corr}";
759  nbins[dim] = nPtBins;
760  min[dim] = 0;
761  max[dim] = fMaxPt;
762  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
763  dim++;
764 
765  title[dim] = "Particle type";
766  nbins[dim] = nParticleTypes;
767  min[dim] = -0.5;
768  max[dim] = 7.5;
769  binEdges[dim] = particleTypeBins;
770  dim++;
771 
772  title[dim] = "N";
773  nbins[dim] = 30;
774  min[dim] = -0.5;
775  max[dim] = 29.5;
776  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
777  dim++;
778 
779  title[dim] = "#it{p}_{T,sum} (GeV)";
780  nbins[dim] = fNPtHistBins;
781  binEdges[dim] = fPtHistBins;
782  min[dim] = fPtHistBins[0];
783  max[dim] = fPtHistBins[fNPtHistBins];
784  dim++;
785 
786  thnname = "JetPerformance/hJetComposition";
787  THnSparse* thn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
788  for (Int_t i = 0; i < dim; i++) {
789  thn->GetAxis(i)->SetTitle(title[i]);
790  thn->SetBinEdges(i, binEdges[i]);
791  }
792 
793  // Hadronic calo energy in each jet
794 
795  // Jet pT vs. Summed energy of hadronic clusters without a matched track
796  histname = "JetPerformance/hHadCaloEnergyUnmatched";
797  htitle = histname + ";#it{p}_{T,jet} (GeV);#it{p}_{T,had} (GeV)";
798  fHistManager.CreateTH2(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins);
799 
800  // Jet pT vs. Summed energy of hadronic clusters with a matched track (before hadronic correction)
801  histname = "JetPerformance/hHadCaloEnergyMatchedNonlincorr";
802  htitle = histname + ";#it{p}_{T,jet} (GeV);#it{p}_{T,had} (GeV)";
803  fHistManager.CreateTH2(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins);
804 
805  // Jet pT vs. Summed energy of hadronic clusters with a matched track (after hadronic correction)
806  histname = "JetPerformance/hHadCaloEnergyMatchedHadCorr";
807  htitle = histname + ";#it{p}_{T,jet} (GeV);#it{p}_{T,had} (GeV)";
808  fHistManager.CreateTH2(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins);
809 
810 }
811 
817 {
819 
820  fNeedEmcalGeom = kTRUE;
821 
822  // Load the PHOS geometry
823  fPHOSGeo = AliPHOSGeometry::GetInstance();
824  if (fPHOSGeo) {
825  AliInfo("Found instance of PHOS geometry!");
826  }
827  else {
828  AliInfo("Creating PHOS geometry!");
829  Int_t runNum = InputEvent()->GetRunNumber();
830  if(runNum<209122) //Run1
831  fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP");
832  else
833  fPHOSGeo = AliPHOSGeometry::GetInstance("Run2");
834 
835  if (fPHOSGeo) {
836  AliOADBContainer geomContainer("phosGeo");
837  geomContainer.InitFromFile("$ALICE_PHYSICS/OADB/PHOS/PHOSMCGeometry.root","PHOSMCRotationMatrixes");
838  TObjArray* matrixes = (TObjArray*)geomContainer.GetObject(runNum,"PHOSRotationMatrixes");
839  for(Int_t mod=0; mod<6; mod++) {
840  if(!matrixes->At(mod)) continue;
841  fPHOSGeo->SetMisalMatrix(((TGeoHMatrix*)matrixes->At(mod)),mod);
842  printf(".........Adding Matrix(%d), geo=%p\n",mod,fPHOSGeo);
843  ((TGeoHMatrix*)matrixes->At(mod))->Print();
844  }
845  }
846  }
847 }
848 
853 {
854  if (fUseAliEventCuts) {
855  if (!fEventCuts.AcceptEvent(InputEvent()))
856  {
857  PostData(1, fOutput);
858  return kFALSE;
859  }
860  }
861  else {
863  }
864  return kTRUE;
865 }
866 
875 {
876  return kTRUE;
877 }
878 
886 {
887 
890  }
891 
892  if (fPlotNeutralJets) {
894  }
895 
896  if (fPlotClustersInJets) {
898  }
899 
900  if (fPlotCellHistograms) {
902  }
903 
904  if (fPlotJetPerformance) {
906  }
907 
908  return kTRUE;
909 }
910 
911 /*
912  * This function fills the histograms for the calorimeter performance study.
913  */
915 {
916  // Define some vars
917  TString histname;
918  Double_t Enonlin;
919  Double_t Ehadcorr;
920  Int_t absId;
921  Double_t ecell;
922  Double_t leadEcell;
923  Int_t leadAbsId;
924 
925  // Get cells from event
926  fCaloCells = InputEvent()->GetEMCALCells();
927  AliVCaloCells* phosCaloCells = InputEvent()->GetPHOSCells();
928 
929  // If MC, get the MC event
930  const AliMCEvent* mcevent = nullptr;
931  if (fGeneratorLevel) {
932  mcevent = MCEvent();
933  }
934 
935  // Loop through clusters and plot cluster THnSparse (centrality, cluster type, E, E-hadcorr, has matched track, M02, Ncells)
936  AliClusterContainer* clusters = 0;
937  const AliVCluster* clus;
938  TString clustersName;
939  TIter nextClusColl(&fClusterCollArray);
940  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
942  clustersName = clusters->GetArrayName();
943  for (AliClusterIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
944 
945  clus = it->second;
946 
947  // Determine cluster type (EMCal/DCal/Phos)
948  ClusterType clusType = kNA;
949  if (clus->IsEMCAL()) {
950  Double_t phi = it->first.Phi_0_2pi();
951  Int_t isDcal = Int_t(phi > fgkEMCalDCalPhiDivide);
952  if (isDcal == 0) {
953  clusType = kEMCal;
954  } else if (isDcal == 1) {
955  clusType = kDCal;
956  }
957  } else if (clus->GetType() == AliVCluster::kPHOSNeutral){
958  clusType = kPHOS;
959  }
960 
961  // rejection reason plots, to make efficiency correction
962  if (clus->IsEMCAL()) {
963  histname = TString::Format("%s/hClusterRejectionReasonEMCal", clusters->GetArrayName().Data());
964  UInt_t rejectionReason = 0;
965  if (!clusters->AcceptCluster(it.current_index(), rejectionReason)) {
966  fHistManager.FillTH2(histname, clusters->GetRejectionReasonBitPosition(rejectionReason), it->first.E());
967  continue;
968  }
969  } else if (clus->GetType() == AliVCluster::kPHOSNeutral){
970  histname = TString::Format("%s/hClusterRejectionReasonPHOS", clusters->GetArrayName().Data());
971  UInt_t rejectionReason = 0;
972  if (!clusters->AcceptCluster(it.current_index(), rejectionReason)) {
973  fHistManager.FillTH2(histname, clusters->GetRejectionReasonBitPosition(rejectionReason), it->first.E());
974  continue;
975  }
976  } else {
977  continue; // avoid CPV clusters
978  }
979 
980  // Fill cluster spectra by SM, and fill cell histograms
981  Enonlin = 0;
982  Ehadcorr = 0;
983  if (clus->IsEMCAL()) {
984 
985  Ehadcorr = clus->GetHadCorrEnergy();
986  Enonlin = clus->GetNonLinCorrEnergy();
988  Enonlin = clus->E();
989  }
990 
991  if (fPlotExotics) {
992  histname = TString::Format("%s/hFcrossEMCal", clusters->GetArrayName().Data());
993  Double_t Fcross = GetFcross(clus, fCaloCells);
994  fHistManager.FillTH3(histname, fCent, Fcross, clus->E());
995  }
996 
997  Int_t sm = fGeom->GetSuperModuleNumber(clus->GetCellAbsId(0));
998  if (sm >=0 && sm < 20) {
999  histname = TString::Format("%s/BySM/hEmcalClusEnergy_SM%d", clusters->GetArrayName().Data(), sm);
1000  fHistManager.FillTH1(histname, clus->E());
1001  }
1002  else {
1003  AliError(Form("Supermodule %d does not exist!", sm));
1004  }
1005 
1006  // Get cells from each accepted cluster, and plot centrality vs. cell energy vs. cell type
1007  histname = TString::Format("Cells/hCellEnergyAccepted");
1008  leadEcell = 0;
1009  for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++){
1010  absId = clus->GetCellAbsId(iCell);
1011  ecell = fCaloCells->GetCellAmplitude(absId);
1012  fHistManager.FillTH3(histname, ecell, fCent, kEMCal); // Note: I don't distinguish EMCal from DCal cells
1013  if (ecell > leadEcell) {
1014  leadEcell = ecell;
1015  leadAbsId = absId;
1016  }
1017  }
1018  // Plot also the leading cell
1019  histname = TString::Format("Cells/hCellEnergyLeading");
1020  fHistManager.FillTH3(histname, leadEcell, fCent, kEMCal);
1021 
1022  } else if (clus->GetType() == AliVCluster::kPHOSNeutral){
1023 
1024  Ehadcorr = clus->GetCoreEnergy();
1025  Enonlin = clus->E();
1026 
1027  Int_t relid[4];
1028  if (fPHOSGeo) {
1029  fPHOSGeo->AbsToRelNumbering(clus->GetCellAbsId(0), relid);
1030  Int_t sm = relid[0];
1031  if (sm >=1 && sm < 5) {
1032  histname = TString::Format("%s/BySM/hPhosClusEnergy_SM%d", clusters->GetArrayName().Data(), sm);
1033  fHistManager.FillTH1(histname, clus->E());
1034  }
1035  else {
1036  AliError(Form("Supermodule %d does not exist!", sm));
1037  }
1038  }
1039 
1040  // Get cells from each accepted cluster, and plot centrality vs. cell energy vs. cell type
1041  histname = TString::Format("Cells/hCellEnergyAccepted");
1042  leadEcell = 0;
1043  for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++){
1044  absId = clus->GetCellAbsId(iCell);
1045  ecell = phosCaloCells->GetCellAmplitude(absId);
1046  fHistManager.FillTH3(histname, ecell, fCent, kPHOS);
1047  if (ecell > leadEcell) {
1048  leadEcell = ecell;
1049  }
1050  }
1051  // Plot also the leading cell
1052  histname = TString::Format("Cells/hCellEnergyLeading");
1053  fHistManager.FillTH3(histname, leadEcell, fCent, kPHOS);
1054  }
1055 
1056  // Check if the cluster has a matched track
1057  Int_t hasMatchedTrack = -1;
1058  Int_t nMatchedTracks = clus->GetNTracksMatched();
1059  if (nMatchedTracks == 0) {
1060  hasMatchedTrack = 0;
1061  } else if (nMatchedTracks > 0) {
1062  hasMatchedTrack = 1;
1063  }
1064 
1065  // Check if the cluster passes the dispersion cut for photon-like cluster (meaningful only for PHOS)
1066  Int_t passedDispersionCut = 0;
1067  if (clus->Chi2() < 2.5*2.5) {
1068  passedDispersionCut = 1;
1069  }
1070 
1071  // Fill info about the cluster
1072  Double_t eta = it->first.Eta();
1073  Double_t phi = it->first.Phi_0_2pi();
1074  Double_t M02 = clus->GetM02();
1075  Int_t nCells = clus->GetNCells();
1076  Double_t distNN = FindNearestNeighborDistance(it->first);
1077 
1078  // If cluster is EMCal, find whether the eta column is even or odd
1079  Int_t isOddEta = -1;
1080  Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta;
1081  if (clus->IsEMCAL()) {
1082  fGeom->GetCellIndex(leadAbsId, nSupMod, nModule, nIphi, nIeta);
1083  fGeom->GetCellPhiEtaIndexInSModule(nSupMod, nModule, nIphi, nIeta, iphi, ieta);
1084  isOddEta = ieta % 2;
1085  }
1086 
1087  // If MC, determine the particle type
1088  // Each detector-level cluster contains an array of labels of each truth-level particle contributing to the cluster
1089  ParticleType particleType1 = kUndefined;
1090  ParticleType particleType2 = kUndefined;
1091 
1092  if (fGeneratorLevel) {
1093 
1094  Int_t label = TMath::Abs(clus->GetLabel()); // returns mc-label of particle that deposited the most energy in the cluster
1095  if (label > 0) { // if the particle has a truth-level match, the label is > 0
1096 
1097  // Method 1: Use AliMCAnalysisUtils to identify all particles
1098  particleType1 = GetParticleType1(clus, mcevent, clusters->GetArray());
1099 
1100  // Method 2: Use MC particle container (and only AliMCAnalysisUtils to find merged pi0)
1101  particleType2 = GetParticleType2(clus, mcevent, label, clusters);
1102 
1103  }
1104  }
1105 
1106  // Standard option: fill once per cluster
1108  FillClusterTHnSparse(clustersName, eta, phi, Enonlin, Ehadcorr, hasMatchedTrack, M02, nCells, passedDispersionCut, distNN, isOddEta, particleType1, particleType2);
1109  }
1110 
1111  if (fPlotCaloCentrality) {
1112 
1113  Double_t eCellCone = 0.;
1114  Int_t nCellsCone = 0;
1115  Double_t eCellSM = 0.;
1116  Int_t nCellsSM = 0;
1117  Double_t areaCone = TMath::Pi() * 0.07 * 0.07;
1118  Double_t areaCell;
1119  Double_t eDensityCone;
1120 
1121  // Get the SM number
1122  Int_t sm = -1;
1123  if (clusType == kEMCal) {
1124  sm = fGeom->GetSuperModuleNumber(clus->GetCellAbsId(0));
1125  areaCell = 0.014*0.014;
1126  }
1127  if (clusType == kPHOS) {
1128  Int_t relid[4];
1129  fPHOSGeo->AbsToRelNumbering(clus->GetCellAbsId(0), relid);
1130  sm = relid[0];
1131  areaCell = 0.014*0.014*(2.2/6.0)*(2.2/6.0); // approximating same r
1132  }
1133 
1134  // Only fill the THnSparse if the cluster is located in a full SM of EMCal or PHOS
1135  if ( (clusType == kEMCal && sm < 10 ) || (clusType == kPHOS && sm < 4) ) {
1136 
1137  eCellCone = GetConeCellEnergy(eta, phi, 0.07) - Enonlin;
1138  nCellsCone = (Int_t)GetConeCellEnergy(eta, phi, 0.07, kTRUE) - nCells;
1139 
1140  eDensityCone = eCellCone / (areaCone - nCells*areaCell);
1141 
1142  if (fPlotCellSMDensity) {
1143  eCellSM = GetSMCellEnergy(sm, clusType) - Enonlin;
1144  nCellsSM = (Int_t)GetSMCellEnergy(sm, clusType, kTRUE) - nCells;
1145  }
1146 
1147  FillClusterTHnSparse(clustersName, eta, phi, Enonlin, eDensityCone, eCellSM, nCellsCone, nCellsSM);
1148 
1149  }
1150 
1151  }
1152 
1153  // If cluster cone option enabled, fill for each R and cone type
1154  if (fPlotClusterCone) {
1155 
1156  // cluster cone, R=0.05
1157  FillClusterTHnSparse(clustersName, eta, phi, Enonlin, Ehadcorr, hasMatchedTrack, M02, nCells, passedDispersionCut, distNN, isOddEta, -1, -1, 0, 0.05, GetConeClusterEnergy(eta, phi, 0.05));
1158  // cluster cone, R=0.1
1159  FillClusterTHnSparse(clustersName, eta, phi, Enonlin, Ehadcorr, hasMatchedTrack, M02, nCells, passedDispersionCut, distNN, isOddEta, -1, -1, 0, 0.1, GetConeClusterEnergy(eta, phi, 0.1));
1160  // cell cone, R=0.05
1161  FillClusterTHnSparse(clustersName, eta, phi, Enonlin, Ehadcorr, hasMatchedTrack, M02, nCells, passedDispersionCut, distNN, isOddEta, -1, -1, 1, 0.05, GetConeCellEnergy(eta, phi, 0.05));
1162  // cell cone, R=0.1
1163  FillClusterTHnSparse(clustersName, eta, phi, Enonlin, Ehadcorr, hasMatchedTrack, M02, nCells, passedDispersionCut, distNN, isOddEta, -1, -1, 1, 0.1, GetConeCellEnergy(eta, phi, 0.1));
1164 
1165  }
1166 
1167  }
1168  }
1169 }
1170 
1171 /*
1172  * This function fills the histograms for the jet performance study.
1173  */
1175 {
1176  TString histname;
1177 
1178  // If MC, get the MC event
1179  const AliMCEvent* mcevent = nullptr;
1180  if (fGeneratorLevel) {
1181  mcevent = MCEvent();
1182  }
1183  else {
1184  return;
1185  }
1186 
1187  // Loop through clusters, and plot M02 for each particle type
1190  const AliVCluster* clus;
1191  for (AliClusterIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
1192 
1193  clus = it->second;
1194 
1195  // Include only EMCal clusters (reject DCal and PHOS clusters)
1196  if (!clus->IsEMCAL()) {
1197  continue;
1198  }
1199  if (it->first.Phi_0_2pi() > fgkEMCalDCalPhiDivide) {
1200  continue;
1201  }
1202 
1203  // If MC, determine the particle type
1204  // Each detector-level cluster contains an array of labels of each truth-level particle contributing to the cluster
1205  ParticleType particleType1 = kUndefined;
1206 
1207  Int_t label = TMath::Abs(clus->GetLabel()); // returns mc-label of particle that deposited the most energy in the cluster
1208  if (label > 0) { // if the particle has a truth-level match, the label is > 0
1209 
1210  // Method 1: Use AliMCAnalysisUtils to identify all particles
1211  particleType1 = GetParticleType1(clus, mcevent, clusters->GetArray());
1212 
1213  }
1214 
1215  // (M02, Eclus, part type)
1216  histname = "JetPerformance/hM02VsParticleType";
1217  fHistManager.FillTH3(histname, clus->GetM02(), clus->GetNonLinCorrEnergy(), particleType1);
1218 
1219  }
1220 
1221  // Loop through jets, to fill various histograms
1222  AliJetContainer* jets = GetJetContainer(0); // there is only a single, det-level jet finder here
1223  for (const auto jet : jets->accepted()) {
1224 
1225  Double_t jetPt = GetJetPt(jet, 0);
1226 
1227  // Keep track of hadronic calo energy in each jet
1228  Double_t hadCaloEnergyUnmatched = 0;
1229  Double_t hadCaloEnergyMatchedNonlincorr = 0;
1230  Double_t hadCaloEnergyMatchedHadCorr = 0;
1231 
1232  // Loop through clusters in each jet, to plot several histograms
1233  Int_t nClusters = jet->GetNumberOfClusters();
1234  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
1235 
1236  clus = jet->Cluster(iClus);
1237 
1238  // Get the particle type of the cluster
1239  ParticleType particleType1 = kUndefined;
1240  Int_t label = TMath::Abs(clus->GetLabel());
1241  if (label > 0) {
1242  particleType1 = GetParticleType1(clus, mcevent, clusters->GetArray());
1243  }
1244 
1245  // Plot M02 for each particle type
1246  histname = "JetPerformance/hM02VsParticleTypeJets";
1247  Double_t x[4] = {clus->GetM02(), clus->GetNonLinCorrEnergy(), particleType1, jetPt};
1248  fHistManager.FillTHnSparse(histname, x);
1249 
1250  // If the cluster is a hadron, sum its energy to compute the jet's hadronic calo energy
1251  if (particleType1 == kHadron) {
1252  Bool_t hasMatchedTrack = (clus->GetNTracksMatched() > 0);
1253  if (hasMatchedTrack) {
1254  hadCaloEnergyMatchedNonlincorr += clus->GetNonLinCorrEnergy();
1255  hadCaloEnergyMatchedHadCorr += clus->GetHadCorrEnergy();
1256  }
1257  else {
1258  hadCaloEnergyUnmatched += clus->GetNonLinCorrEnergy();
1259  }
1260  }
1261 
1262  }
1263 
1264  // Fill hadronic calo energy in each jet
1265 
1266  // (Jet pT, Summed energy of hadronic clusters without a matched track)
1267  histname = "JetPerformance/hHadCaloEnergyUnmatched";
1268  fHistManager.FillTH2(histname, jetPt, hadCaloEnergyUnmatched);
1269 
1270  // (Jet pT vs. Summed energy of hadronic clusters with a matched track (before hadronic correction))
1271  histname = "JetPerformance/hHadCaloEnergyMatchedNonlincorr";
1272  fHistManager.FillTH2(histname, jetPt, hadCaloEnergyMatchedNonlincorr);
1273 
1274  // (Jet pT vs. Summed energy of hadronic clusters with a matched track (after hadronic correction))
1275  histname = "JetPerformance/hHadCaloEnergyMatchedHadCorr";
1276  fHistManager.FillTH2(histname, jetPt, hadCaloEnergyMatchedHadCorr);
1277 
1278  // Loop through particle types, and plot jet composition for each particle type
1279  histname = "JetPerformance/hJetComposition";
1280  for (Int_t type = 0; type < 8; type++) {
1281 
1282  ParticleType particleType1 = kUndefined;
1283  Double_t nSum = 0;
1284  Double_t pTsum = 0;
1285 
1286  // Loop through clusters in jet, and add to sum if cluster matches particle type
1287  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
1288 
1289  clus = jet->Cluster(iClus);
1290 
1291  Int_t label = TMath::Abs(clus->GetLabel());
1292  if (label > 0) {
1293  particleType1 = GetParticleType1(clus, mcevent, clusters->GetArray());
1294  }
1295 
1296  if (type == particleType1) {
1297  nSum++;
1298  pTsum += clus->GetNonLinCorrEnergy();
1299  }
1300  }
1301 
1302  // Fill jet composition histogram
1303  Double_t x[4] = {jetPt, 1.*type, nSum, pTsum};
1304  fHistManager.FillTHnSparse(histname, x);
1305 
1306  }
1307  }
1308 }
1309 
1310 /*
1311  * Compute the MC particle type using AliMCAnalysisUtils
1312  */
1313 AliAnalysisTaskEmcalVsPhos::ParticleType AliAnalysisTaskEmcalVsPhos::GetParticleType1(const AliVCluster* clus, const AliMCEvent* mcevent, const TClonesArray* clusArray)
1314 {
1315  ParticleType particleType = kUndefined;
1316 
1317  AliMCAnalysisUtils mcUtils;
1318  Int_t tag = mcUtils.CheckOrigin(clus->GetLabels(), clus->GetNLabels(), mcevent, clusArray);
1319 
1320  Bool_t isPhoton = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton);
1321  Bool_t isPi0 = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0);
1322  Bool_t isConversion = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion);
1323  Bool_t isEta = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCEta);
1324  Bool_t isPion = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCPion);
1325  Bool_t isKaon = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCKaon);
1326  Bool_t isProton = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCProton);
1327  Bool_t isAntiProton = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCAntiProton);
1328  Bool_t isNeutron = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCNeutron);
1329  Bool_t isAntiNeutron = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCAntiNeutron);
1330  Bool_t isElectron = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron);
1331  Bool_t isMuon = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCMuon);
1332 
1333  if (isPi0) {
1334  if (isConversion) {
1335  particleType = kPi0Conversion;
1336  }
1337  else {
1338  particleType = kPi0;
1339  }
1340  }
1341  else if (isEta) {
1342  particleType = kEta;
1343  }
1344  else if (isPhoton) {
1345  particleType = kPhoton;
1346  }
1347  else if (isPion || isKaon || isProton || isAntiProton || isNeutron || isAntiNeutron) {
1348  particleType = kHadron;
1349  }
1350  else if (isElectron) {
1351  particleType = kElectron;
1352  }
1353  else if (isMuon) {
1354  particleType = kMuon;
1355  }
1356  else {
1357  particleType = kOther;
1358  }
1359  return particleType;
1360 }
1361 
1362 /*
1363  * Compute the MC particle type using the MC particle container (and only AliMCAnalysisUtils to find merged pi0)
1364  */
1365 AliAnalysisTaskEmcalVsPhos::ParticleType AliAnalysisTaskEmcalVsPhos::GetParticleType2(const AliVCluster* clus, const AliMCEvent* mcevent, Int_t label, const AliClusterContainer* clusters)
1366 {
1367  ParticleType particleType = kUndefined;
1368 
1369  AliAODMCParticle *part = fGeneratorLevel->GetMCParticleWithLabel(label);
1370  if (part) {
1371 
1372  TString histname = TString::Format("%s/hClusterRejectionReasonMC", clusters->GetArrayName().Data());
1373  UInt_t rejectionReason = 0;
1374  if (!fGeneratorLevel->AcceptMCParticle(part, rejectionReason)) {
1375  fHistManager.FillTH2(histname, fGeneratorLevel->GetRejectionReasonBitPosition(rejectionReason), clus->E());
1376  return particleType;
1377  }
1378 
1379  if (part->GetGeneratorIndex() == 0) { // generator index in cocktail
1380 
1381  // select charged pions, protons, kaons, electrons, muons
1382  Int_t pdg = TMath::Abs(part->PdgCode()); // abs value ensures both particles and antiparticles are included
1383 
1384  if (pdg == 22) { // gamma 22
1385 
1386  AliMCAnalysisUtils mcUtils;
1387  Int_t tag;
1388  mcUtils.CheckOverlapped2GammaDecay(clus->GetLabels(), clus->GetNLabels(), part->GetMother(), mcevent, tag);
1389 
1390  if (mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0)) {
1391  if (mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)) {
1392  particleType = kPi0Conversion;
1393  }
1394  else {
1395  particleType = kPi0;
1396  }
1397  }
1398  else if (mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)) {
1399  particleType = kEta;
1400  }
1401  else { // direct photon
1402  particleType = kPhoton;
1403  }
1404 
1405  }
1406  else if (pdg == 211 || 2212 || 321 || 2112) { // pi+ 211, proton 2212, K+ 321, neutron 2112
1407  particleType = kHadron;
1408  }
1409  else if (pdg == 11) { // e- 11
1410  particleType = kElectron;
1411  }
1412  else if (pdg == 13) { // mu- 13
1413  particleType = kMuon;
1414  }
1415  else {
1416  particleType = kOther;
1417  }
1418  }
1419  }
1420  return particleType;
1421 }
1422 
1423 
1424 /*
1425  * This function fills the cluster THnSparse.
1426  */
1427 void AliAnalysisTaskEmcalVsPhos::FillClusterTHnSparse(TString clustersName, Double_t eta, Double_t phi, Double_t Enonlin, Double_t Ehadcorr, Int_t hasMatchedTrack, Double_t M02, Int_t nCells, Int_t passedDispersionCut, Double_t distNN, Int_t isOddEta, Int_t particleType1, Int_t particleType2, Int_t coneType, Double_t R, Double_t Econe)
1428 {
1429  Double_t contents[30]={0};
1430  TString histname = TString::Format("%s/clusterObservables", clustersName.Data());
1431  THnSparse* histClusterObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1432  if (!histClusterObservables) return;
1433  for (Int_t i = 0; i < histClusterObservables->GetNdimensions(); i++) {
1434  TString title(histClusterObservables->GetAxis(i)->GetTitle());
1435  if (title=="Centrality %")
1436  contents[i] = fCent;
1437  else if (title=="#eta")
1438  contents[i] = eta;
1439  else if (title=="#phi")
1440  contents[i] = phi;
1441  else if (title=="#it{E}_{clus} (GeV)")
1442  contents[i] = Enonlin;
1443  else if (title=="#it{E}_{clus, hadcorr} or #it{E}_{core} (GeV)")
1444  contents[i] = Ehadcorr;
1445  else if (title=="Matched track")
1446  contents[i] = hasMatchedTrack;
1447  else if (title=="M02")
1448  contents[i] = M02;
1449  else if (title=="Ncells")
1450  contents[i] = nCells;
1451  else if (title=="Dispersion cut")
1452  contents[i] = passedDispersionCut;
1453  else if (title=="#DeltaR_{NN}")
1454  contents[i] = distNN;
1455  else if (title=="Even/odd eta")
1456  contents[i] = isOddEta;
1457  else if (title=="Particle type1")
1458  contents[i] = particleType1;
1459  else if (title=="Particle type2")
1460  contents[i] = particleType2;
1461  else if (title=="Cone type")
1462  contents[i] = coneType;
1463  else if (title=="R")
1464  contents[i] = R;
1465  else if (title=="#it{E}_{cone} (GeV)")
1466  contents[i] = Econe;
1467  else
1468  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1469  }
1470  histClusterObservables->Fill(contents);
1471 
1472 }
1473 
1474 /*
1475  * This function fills the cluster THnSparse (alternate signature, used for local density option).
1476  */
1477 void AliAnalysisTaskEmcalVsPhos::FillClusterTHnSparse(TString clustersName, Double_t eta, Double_t phi, Double_t Enonlin, Double_t eCellCone, Double_t eCellSM, Int_t nCellsCone, Int_t nCellsSM)
1478 {
1479  Double_t contents[30]={0};
1480  TString histname = TString::Format("%s/clusterObservables", clustersName.Data());
1481  THnSparse* histClusterObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1482  if (!histClusterObservables) return;
1483  for (Int_t i = 0; i < histClusterObservables->GetNdimensions(); i++) {
1484  TString title(histClusterObservables->GetAxis(i)->GetTitle());
1485  if (title=="Centrality %")
1486  contents[i] = fCent;
1487  else if (title=="#eta")
1488  contents[i] = eta;
1489  else if (title=="#phi")
1490  contents[i] = phi;
1491  else if (title=="#it{E}_{clus} (GeV)")
1492  contents[i] = Enonlin;
1493  else if (title=="#it{#rho}_{cell cone} (GeV)")
1494  contents[i] = eCellCone;
1495  else if (title=="Ncells cone")
1496  contents[i] = nCellsCone;
1497  else if (title=="#it{E}_{cell SM} (GeV)")
1498  contents[i] = eCellSM;
1499  else if (title=="Ncells SM")
1500  contents[i] = nCellsSM;
1501  else
1502  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1503  }
1504  histClusterObservables->Fill(contents);
1505 
1506 }
1507 
1508 /*
1509  * This function fills cell histograms.
1510  */
1512 {
1513  // Define some vars
1514  TString histname;
1515  Int_t absId;
1516  Double_t ecell;
1517 
1518  // Get cells from event
1519  fCaloCells = InputEvent()->GetEMCALCells();
1520  AliVCaloCells* phosCaloCells = InputEvent()->GetPHOSCells();
1521 
1522  // Loop through all cells and fill histos
1523  Int_t sm;
1524  Int_t relid[4];
1525  Double_t patchSumEMCal[20] = {0.};
1526  Double_t patchSumPHOS[4] = {0.};
1527  for (Int_t i=0; i<fCaloCells->GetNumberOfCells(); i++) {
1528 
1529  absId = fCaloCells->GetCellNumber(i);
1530  ecell = fCaloCells->GetCellAmplitude(absId);
1531 
1532  // Fill cell histo
1533  histname = TString::Format("Cells/hCellEnergyAll");
1534  fHistManager.FillTH3(histname, ecell, fCent, kEMCal); // Note: I don't distinguish EMCal from DCal cells
1535 
1536  // Fill cell patch histo, per SM
1537  sm = fGeom->GetSuperModuleNumber(absId);
1538  if (sm >=0 && sm < 20) {
1539  patchSumEMCal[sm] += ecell;
1540  }
1541 
1542  }
1543 
1544  for (Int_t i=0; i<phosCaloCells->GetNumberOfCells(); i++) {
1545 
1546  absId = phosCaloCells->GetCellNumber(i);
1547  ecell = phosCaloCells->GetCellAmplitude(absId);
1548 
1549  if (absId < 0) {
1550  continue; // skip CPV cells
1551  }
1552 
1553  // Fill cell histo
1554  histname = TString::Format("Cells/hCellEnergyAll");
1555  fHistManager.FillTH3(histname, ecell, fCent, kPHOS);
1556 
1557  // Fill cell patch histo, per SM
1558  fPHOSGeo->AbsToRelNumbering(absId, relid);
1559  sm = relid[0];
1560  if (sm >=1 && sm < 5) {
1561  patchSumPHOS[sm-1] += ecell;
1562  }
1563 
1564  }
1565 
1566  for (Int_t sm = 0; sm < 20; sm++) {
1567  histname = TString::Format("Cells/BySM/hEmcalPatchEnergy_SM%d", sm);
1568  fHistManager.FillTH2(histname, patchSumEMCal[sm], fCent);
1569  }
1570 
1571  for (Int_t sm = 1; sm < 5; sm++) {
1572  histname = TString::Format("Cells/BySM/hPhosPatchEnergy_SM%d", sm);
1573  fHistManager.FillTH2(histname, patchSumPHOS[sm-1], fCent);
1574  }
1575 
1576 }
1577 
1578 /*
1579  * This function fills neutral jet histograms.
1580  */
1582 {
1583  AliJetContainer* jets = 0;
1584  TIter nextJetColl(&fJetCollArray);
1585  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
1586 
1587  // plot neutral jets THnSparse (centrality, eta, phi, E, Nclusters)
1588  Double_t contents[30]={0};
1589  TString histname = TString::Format("%s/hClusterJetObservables", jets->GetArrayName().Data());
1590  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1591  if (!histJetObservables) return;
1592 
1593  for (auto jet : jets->accepted()) {
1594 
1595  for (Int_t i = 0; i < histJetObservables->GetNdimensions(); i++) {
1596  TString title(histJetObservables->GetAxis(i)->GetTitle());
1597  if (title=="Centrality (%)")
1598  contents[i] = fCent;
1599  else if (title=="#eta_{jet}")
1600  contents[i] = jet->Eta();
1601  else if (title=="#phi_{jet} (rad)")
1602  contents[i] = jet->Phi_0_2pi();
1603  else if (title=="#it{E}_{T} (GeV)")
1604  contents[i] = jet->Pt();
1605  else if (title=="#rho (GeV/#it{c})")
1606  contents[i] = jet->Pt() / jet->Area();
1607  else if (title=="N_{clusters}")
1608  contents[i] = jet->GetNumberOfClusters();
1609  else if (title=="#it{E}_{T}, acc clus within R (GeV)")
1610  contents[i] = GetConeClusterEnergy(jet->Eta(), jet->Phi_0_2pi(), jets->GetJetRadius());
1611  else if (title=="#it{E}_{T}, acc cell within R (GeV)")
1612  contents[i] = GetConeCellEnergy(jet->Eta(), jet->Phi_0_2pi(), jets->GetJetRadius());
1613  else
1614  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1615  }
1616  histJetObservables->Fill(contents);
1617 
1618  }
1619  }
1620 }
1621 
1622 /*
1623  * This function fills clusters within jets and estimates the JES shift due to the bump.
1624  */
1626 {
1627  TString histname;
1628  Double_t rho;
1629  AliJetContainer* jets = 0;
1630  TIter nextJetColl(&fJetCollArray);
1631  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
1632 
1633  rho = jets->GetRhoVal();
1634 
1635  for (const auto jet : jets->accepted()) {
1636 
1637  // Fill cluster spectra of clusters within jets
1638  //(centrality, cluster energy, jet pT, jet eta, jet phi)
1639  histname = TString::Format("%s/hClustersInJets", jets->GetArrayName().Data());
1640  Int_t nClusters = jet->GetNumberOfClusters();
1641  AliVCluster* clus;
1642  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
1643 
1644  clus = jet->Cluster(iClus);
1645 
1647  Double_t x[5] = {fCent, clus->E(), GetJetPt(jet, rho), jet->Eta(), jet->Phi_0_2pi()};
1648  fHistManager.FillTHnSparse(histname, x);
1649  }
1650  else {
1651  Double_t x[4] = {clus->E(), GetJetPt(jet, rho), jet->Eta(), jet->Phi_0_2pi()};
1652  fHistManager.FillTHnSparse(histname, x);
1653  }
1654 
1655  }
1656 
1657  // Loop through clusters, and plot estimated shift in JES due to cluster bump
1658  // Only do for 0-10% centrality, and for EMCal/DCal
1659  Double_t eclus;
1660  Double_t shift;
1661  Double_t shiftSum = 0;
1663  if (GetJetType(jet) > -0.5 && GetJetType(jet) < 1.5) {
1664  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
1665  clus = jet->Cluster(iClus);
1666  eclus = clus->E();
1667  if (eclus > 0.5) {
1668  shift = 0.79 * TMath::Exp(-0.5 * ((eclus - 3.81) / 1.50)*((eclus - 3.81) / 1.50) );
1669  shiftSum += shift;
1670  }
1671  }
1672  histname = TString::Format("%s/hCaloJESshift", jets->GetArrayName().Data());
1673  fHistManager.FillTH3(histname, GetJetType(jet), GetJetPt(jet, rho), shiftSum);
1674  }
1675  }
1676 
1677  }
1678  }
1679 }
1680 
1685 {
1686  Double_t pT = jet->Pt() - rho * jet->Area();
1687  return pT;
1688 }
1689 
1694 {
1695  Double_t deltaPhi = TMath::Abs(part.Phi_0_2pi() - phiRef);
1696  Double_t deltaEta = TMath::Abs(part.Eta() - etaRef);
1697  Double_t deltaR = TMath::Sqrt( deltaPhi*deltaPhi + deltaEta*deltaEta );
1698  return deltaR;
1699 }
1700 
1705 {
1706  Double_t deltaPhi = TMath::Abs(phi1-phi2);
1707  Double_t deltaEta = TMath::Abs(eta1-eta2);
1708  Double_t deltaR = TMath::Sqrt( deltaPhi*deltaPhi + deltaEta*deltaEta );
1709  return deltaR;
1710 }
1711 
1716 {
1717  UInt_t jetType = jet->GetJetAcceptanceType();
1718  Double_t type = -1;
1719  if (jetType & AliEmcalJet::kEMCAL) {
1720  type = kEMCal;
1721  }
1722  else if (jetType & AliEmcalJet::kDCALonly) {
1723  type = kDCal;
1724  }
1725  else if (jetType & AliEmcalJet::kPHOS) {
1726  type = kPHOS;
1727  }
1728 
1729  return type;
1730 }
1731 
1735 Double_t AliAnalysisTaskEmcalVsPhos::GetFcross(const AliVCluster *cluster, AliVCaloCells *cells)
1736 {
1737  Int_t AbsIdseed = -1;
1738  Double_t Eseed = 0;
1739  for (Int_t i = 0; i < cluster->GetNCells(); i++) {
1740  if (cells->GetCellAmplitude(cluster->GetCellAbsId(i)) > Eseed) {
1741  Eseed = cells->GetCellAmplitude(cluster->GetCellAbsId(i));
1742  AbsIdseed = cluster->GetCellAbsId(i);
1743  }
1744  }
1745 
1746  if (Eseed < 1e-9) {
1747  return 100;
1748  }
1749 
1750  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
1751  fGeom->GetCellIndex(AbsIdseed,imod,iTower,iIphi,iIeta);
1752  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi,iIeta,iphi,ieta);
1753 
1754  //Get close cells index and energy, not in corners
1755 
1756  Int_t absID1 = -1;
1757  Int_t absID2 = -1;
1758 
1759  if (iphi < AliEMCALGeoParams::fgkEMCALRows-1) {
1760  absID1 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
1761  }
1762  if (iphi > 0) {
1763  absID2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
1764  }
1765 
1766  // In case of cell in eta = 0 border, depending on SM shift the cross cell index
1767 
1768  Int_t absID3 = -1;
1769  Int_t absID4 = -1;
1770 
1771  if (ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2)) {
1772  absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
1773  absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
1774  }
1775  else if (ieta == 0 && imod%2) {
1776  absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
1777  absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
1778  }
1779  else {
1780  if (ieta < AliEMCALGeoParams::fgkEMCALCols-1) {
1781  absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
1782  }
1783  if (ieta > 0) {
1784  absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
1785  }
1786  }
1787 
1788  Double_t ecell1 = cells->GetCellAmplitude(absID1);
1789  Double_t ecell2 = cells->GetCellAmplitude(absID2);
1790  Double_t ecell3 = cells->GetCellAmplitude(absID3);
1791  Double_t ecell4 = cells->GetCellAmplitude(absID4);
1792 
1793  Double_t Ecross = ecell1 + ecell2 + ecell3 + ecell4;
1794 
1795  Double_t Fcross = 1 - Ecross/Eseed;
1796 
1797  return Fcross;
1798 }
1799 
1804 {
1805  Double_t distNN = 10.;
1806  Double_t etaRef = clusterRef.Eta();
1807  Double_t phiRef = clusterRef.Phi_0_2pi();
1808 
1810  AliTLorentzVector clusNNcand;
1811  for (auto clusIterator : clusters->accepted_momentum() ) {
1812 
1813  clusNNcand.Clear();
1814  clusNNcand = clusIterator.first;
1815 
1816  Double_t distNNcand = GetDeltaR(clusNNcand, etaRef, phiRef);
1817 
1818  if (distNNcand < distNN && distNNcand > 0.001) {
1819  distNN = distNNcand;
1820  }
1821 
1822  }
1823 
1824  return distNN;
1825 
1826 }
1827 
1832 {
1834  AliTLorentzVector clus;
1835  Double_t energy = 0.;
1836  for (auto clusIterator : clusCont->accepted_momentum() ) {
1837 
1838  clus.Clear();
1839  clus = clusIterator.first;
1840 
1841  if (GetDeltaR(clus, etaRef, phiRef) < R) {
1842  energy += clus.E();
1843  }
1844  }
1845  return energy;
1846 }
1847 
1853 {
1854  Double_t energy = 0.;
1855  Double_t nCells = 0.;
1856 
1857  // Get cells from event
1858  fCaloCells = InputEvent()->GetEMCALCells();
1859  AliVCaloCells* phosCaloCells = InputEvent()->GetPHOSCells();
1860 
1861  Double_t eta;
1862  Double_t phi;
1863  Int_t absId;
1864  TVector3 pos;
1865  Int_t relid[4];
1866  Int_t sm;
1867 
1868  for (Int_t i=0; i<fCaloCells->GetNumberOfCells(); i++) {
1869 
1870  absId = fCaloCells->GetCellNumber(i);
1871 
1872  fGeom->EtaPhiFromIndex(absId, eta, phi);
1873  phi = TVector2::Phi_0_2pi(phi);
1874 
1875  if (GetDeltaR(eta, phi, etaRef, phiRef) < R) {
1876 
1877  if (fExcludeRejectedCells) {
1878  if (IsCellRejected(absId, kEMCal)) {
1879  continue;
1880  }
1881  }
1882 
1883  energy += fCaloCells->GetCellAmplitude(absId);
1884  nCells += 1.;
1885  }
1886 
1887  }
1888 
1889  for (Int_t i=0; i<phosCaloCells->GetNumberOfCells(); i++) {
1890 
1891  absId = phosCaloCells->GetCellNumber(i);
1892 
1893  fPHOSGeo->AbsToRelNumbering(absId, relid);
1894  sm = relid[0];
1895  if (sm < 1 || sm > 4) {
1896  continue;
1897  }
1898 
1899  fPHOSGeo->RelPosInAlice(absId, pos); // pos then contains (x,y,z) coordinate of cell
1900  eta = pos.Eta();
1901  phi = pos.Phi();
1902  phi = TVector2::Phi_0_2pi(phi);
1903 
1904  if (GetDeltaR(eta, phi, etaRef, phiRef) < R) {
1905 
1906  if (fExcludeRejectedCells) {
1907  if (IsCellRejected(absId, kPHOS)) {
1908  continue;
1909  }
1910  }
1911 
1912  energy += phosCaloCells->GetCellAmplitude(absId);
1913  nCells += 1.;
1914  }
1915 
1916  }
1917 
1918  if (returnNcells) {
1919  return nCells;
1920  }
1921 
1922  return energy;
1923 }
1924 
1930 {
1931  Double_t energy = 0.;
1932  Double_t nCells = 0.;
1933  Int_t absId;
1934  Int_t cellSM;
1935  Int_t relid[4];
1936 
1937  if (clusType == kEMCal) {
1938 
1939  fCaloCells = InputEvent()->GetEMCALCells();
1940 
1941  for (Int_t i=0; i<fCaloCells->GetNumberOfCells(); i++) {
1942 
1943  absId = fCaloCells->GetCellNumber(i);
1944  cellSM = fGeom->GetSuperModuleNumber(absId);
1945 
1946  if (cellSM == sm) {
1947 
1948  if (fExcludeRejectedCells) {
1949  if (IsCellRejected(absId, kEMCal)) {
1950  continue;
1951  }
1952  }
1953 
1954  energy += fCaloCells->GetCellAmplitude(absId);
1955  nCells += 1.;
1956 
1957  }
1958  }
1959  }
1960 
1961  if (clusType == kPHOS) {
1962 
1963  AliVCaloCells* phosCaloCells = InputEvent()->GetPHOSCells();
1964 
1965  for (Int_t i=0; i<phosCaloCells->GetNumberOfCells(); i++) {
1966 
1967  absId = phosCaloCells->GetCellNumber(i);
1968 
1969  fPHOSGeo->AbsToRelNumbering(absId, relid);
1970  cellSM = relid[0];
1971 
1972  if (cellSM == sm) {
1973 
1974  if (fExcludeRejectedCells) {
1975  if (IsCellRejected(absId, kPHOS)) {
1976  continue;
1977  }
1978  }
1979 
1980  energy += phosCaloCells->GetCellAmplitude(absId);
1981  nCells += 1.;
1982 
1983  }
1984  }
1985  }
1986 
1987  if (returnNcells) {
1988  return nCells;
1989  }
1990 
1991  return energy;
1992 }
1993 
1998 {
2000  AliVCluster* cluster;
2001 
2002  UInt_t rejectionReason;
2003  Bool_t skipCell = kFALSE;
2004  Int_t clusType;
2005 
2007  for (AliClusterIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
2008 
2009  if (!clusCont->AcceptCluster(it.current_index(), rejectionReason)) {
2010 
2011  cluster = it->second;
2012 
2013  // check that the cell type and cluster type are the same
2014  if (cluster->IsEMCAL()) {
2015  clusType = kEMCal;
2016  }
2017  if (cluster->GetType() == AliVCluster::kPHOSNeutral) {
2018  clusType = kPHOS;
2019  }
2020  if (clusType != cellType) {
2021  continue;
2022  }
2023 
2024  // skip the cell if it belongs to a rejected cluster
2025  for (Int_t i = 0; i < cluster->GetNCells(); i++) {
2026 
2027  if (absId == cluster->GetCellAbsId(i)) {
2028  skipCell = kTRUE;
2029  }
2030 
2031  }
2032  }
2033  }
2034  return skipCell;
2035 }
2036 
Int_t pdg
void Print(std::ostream &o, const char *name, Double_t dT, Double_t dVM, Double_t alldT, Double_t alldVM)
Definition: PlotSysInfo.C:121
Double_t Area() const
Definition: AliEmcalJet.h:123
TObjArray fClusterCollArray
cluster collection array
Double_t GetRhoVal() const
double Double_t
Definition: External.C:58
ParticleType GetParticleType2(const AliVCluster *clus, const AliMCEvent *mcevent, Int_t label, const AliClusterContainer *clusters)
const char * title
Definition: MakeQAPdf.C:27
AliJetContainer * GetJetContainer(Int_t i=0) const
Bool_t fPlotEvenOddEta
Set whether to add axis to THnSparse separating even/odd eta columns.
PHOS acceptance.
Definition: AliEmcalJet.h:68
Bool_t fPlotClustersInJets
Set whether to plot histogram of clusters within jets.
UInt_t fOffTrigger
offline trigger for event selection
Double_t GetSMCellEnergy(Int_t sm, Int_t clusType, Bool_t returnNcells=kFALSE)
const AliClusterIterableMomentumContainer accepted_momentum() const
bidirectional stl iterator over the EMCAL iterable container
Double_t GetFcross(const AliVCluster *cluster, AliVCaloCells *cells)
Int_t fNCentHistBins
! number of cent bins
energy
Definition: HFPtSpectrum.C:44
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
Double_t * fPtHistBins
! variable pt bins
Declaration of class AliTLorentzVector.
void FillTH3(const char *hname, double x, double y, double z, double weight=1., Option_t *opt="")
Fill a 3D histogram within the container.
Bool_t fPlotCellHistograms
Set whether to plot cell histograms.
Bool_t fUseAliEventCuts
Flag to use AliEventCuts (otherwise AliAnalysisTaskEmcal will be used)
ParticleType GetParticleType1(const AliVCluster *clus, const AliMCEvent *mcevent, const TClonesArray *clusArray)
Int_t fNM02HistBins
! number of variable M02 bins
Double_t FindNearestNeighborDistance(AliTLorentzVector cluster)
Int_t CheckOrigin(Int_t label, const AliMCEvent *mcevent)
THistManager fHistManager
Histogram manager.
AliPHOSGeometry * fPHOSGeo
! phos geometry
virtual Bool_t AcceptMCParticle(const AliAODMCParticle *vp, UInt_t &rejectionReason) const
Pi0 (merged pi0) with conversion of one photon (may be only partially contained in cluster) ...
const Int_t nPtBins
void FillTHnSparse(const char *name, const double *x, double weight=1., Option_t *opt="")
Double_t GetConeCellEnergy(Double_t etaRef, Double_t phiRef, Double_t R, Bool_t returnNcells=kFALSE)
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
Bool_t fPlotStandardClusterTHnSparse
Set whether to plot "standard" axes in cluster THnSparse.
unsigned int UInt_t
Definition: External.C:33
Int_t fNPtHistBins
! number of variable pt bins
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
Bool_t fPlotNeutralJets
Set whether to plot neutral jet histo.
AliEMCALGeometry * fGeom
!emcal geometry
UInt_t GetJetAcceptanceType() const
Definition: AliEmcalJet.h:262
Bool_t IsCellRejected(Int_t absId, Int_t cellType)
Double_t Phi_0_2pi() const
EMCal acceptance.
Definition: AliEmcalJet.h:62
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
Bool_t fPlotClusterCone
Set whether to plot sum of energy surrounding cluster in THnSparse.
Bool_t fPlotFineGrainedEtaPhi
Set whether to plot fine-grained eta-phi bins in cluster THnSparse.
BeamType fForceBeamType
forced beam type
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
virtual Bool_t AcceptCluster(Int_t i, UInt_t &rejectionReason) const
Double_t fCent
!event centrality
Bool_t fExcludeRejectedCells
Set whether to exclude cells from rejected clusters in cone/SM studies.
Bool_t fPlotJetPerformance
Set whether to plot jet reconstruction studies with M02 cut.
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
Double_t GetJetType(const AliEmcalJet *jet)
Bool_t fPlotCellSMDensity
Set whether to plot SM cell density when computing local density.
TObjArray fJetCollArray
jet collection array
Bool_t fPlotClusterHistograms
Set whether to plot cluster histograms.
Bool_t fUseManualEventCuts
Flag to use manual event cuts.
TList * fEventCutList
! Output list for event cut histograms
AliVCaloCells * fCaloCells
!cells
void CheckOverlapped2GammaDecay(const Int_t *labels, Int_t nlabels, Int_t mesonIndex, const AliMCEvent *mcevent, Int_t &tag)
Double_t Pt() const
Definition: AliEmcalJet.h:102
Bool_t fPlotNearestNeighborDistribution
Set whether to plot nearest neighbor axis in cluster THnSparse.
const AliClusterIterableMomentumContainer all_momentum() const
Double_t * fM02HistBins
! variable M02 bins
virtual Bool_t IsEventSelected()
Performing event selection.
static Double_t * GenerateFixedBinArray(Int_t n, Double_t min, Double_t max)
Bool_t fPlotFineGrainedCentrality
Set whether to plot a more fine grained centrality binning.
Float_t GetJetRadius() const
AliEmcalList * fOutput
!output list
AliMCParticleContainer * GetMCParticleContainer(Int_t i=0) const
Double_t GetDeltaR(AliTLorentzVector part, Double_t etaRef, Double_t phiRef)
Double_t GetConeClusterEnergy(Double_t etaRef, Double_t phiRef, Double_t R)
Definition: External.C:220
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
void FillClusterTHnSparse(TString clustersName, Double_t eta, Double_t phi, Double_t Enonlin, Double_t Ehadcorr, Int_t hasMatchedTrack, Double_t M02, Int_t nCells, Int_t passDispersionCut, Double_t distNN, Int_t isOddEta, Int_t particleType1=-1, Int_t particleType2=-1, Int_t coneType=0, Double_t R=0., Double_t Econe=0.)
void SetRejectionReasonLabels(TAxis *axis)
AliMCParticleContainer * fGeneratorLevel
! generator level container
void UserCreateOutputObjects()
Main initialization function on the worker.
const Int_t nbins
const AliJetIterableContainer accepted() const
bool Bool_t
Definition: External.C:53
Double_t GetJetPt(const AliEmcalJet *jet, Double_t rho)
Bool_t fPlotCaloCentrality
Set whether to bin cluster THnSparse in calorimeter local density.
Bool_t fPlotClusWithoutNonLinCorr
If true, use pre-nonlincorr energy in cluster thnsparse.
DCal acceptance – spans ONLY DCal (no PHOS or gap)
Definition: AliEmcalJet.h:66
Class with analysis utils for simulations.
THnSparse * CreateTHnSparse(const char *name, const char *title, int ndim, const int *nbins, const double *min, const double *max, Option_t *opt="")
Create a new THnSparse within the container.
Container structure for EMCAL clusters.
virtual AliAODMCParticle * GetMCParticleWithLabel(Int_t lab) const
Bool_t fNeedEmcalGeom
whether or not the task needs the emcal geometry
Bool_t CheckTagBit(Int_t tag, UInt_t test) const
Bool_t fPlotExotics
Set whether to plot exotic cluster study.
Container for jet within the EMCAL jet framework.
TH3 * CreateTH3(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, int nbinsz, double zmin, double zmax, Option_t *opt="")
Create a new TH2 within the container.
Float_t fMaxPt
Histogram pt limit.
Study of EMCal vs. PHOS clusters.
AliEventCuts fEventCuts
event selection utility
static Double_t fgkEMCalDCalPhiDivide
phi value used to distinguish between DCal and EMCal