AliPhysics  068200c (068200c)
 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 
52 
56 
62  fPlotClusterHistograms(kTRUE),
63  fPlotNeutralJets(kFALSE),
64  fPlotClustersInJets(kFALSE),
65  fPlotCellHistograms(kTRUE),
66  fPlotClusWithoutNonLinCorr(kFALSE),
67  fPlotExotics(kFALSE),
68  fPlotStandardClusterTHnSparse(kTRUE),
69  fPlotNearestNeighborDistribution(kFALSE),
70  fPlotClusterCone(kFALSE),
71  fPlotCaloCentrality(kFALSE),
72  fPlotFineGrainedEtaPhi(kFALSE),
73  fPlotEvenOddEta(kFALSE),
74  fPlotCellSMDensity(kFALSE),
75  fExcludeRejectedCells(kFALSE),
76  fPlotFineGrainedCentrality(kFALSE),
77  fMaxPt(200),
78  fNCentHistBins(0),
79  fCentHistBins(0),
80  fNPtHistBins(0),
81  fPtHistBins(0),
82  fNM02HistBins(0),
83  fM02HistBins(0),
84  fUseAliEventCuts(kTRUE),
85  fEventCuts(0),
86  fEventCutList(0),
87  fUseManualEventCuts(kFALSE),
88  fGeneratorLevel(0),
89  fPHOSGeo(nullptr),
90  fHistManager()
91 {
93 }
94 
101  AliAnalysisTaskEmcalJet(name, kTRUE),
102  fPlotClusterHistograms(kTRUE),
103  fPlotNeutralJets(kFALSE),
104  fPlotClustersInJets(kFALSE),
105  fPlotCellHistograms(kTRUE),
106  fPlotClusWithoutNonLinCorr(kFALSE),
107  fPlotExotics(kFALSE),
108  fPlotStandardClusterTHnSparse(kTRUE),
109  fPlotNearestNeighborDistribution(kFALSE),
110  fPlotClusterCone(kFALSE),
111  fPlotCaloCentrality(kFALSE),
112  fPlotFineGrainedEtaPhi(kFALSE),
113  fPlotEvenOddEta(kFALSE),
114  fPlotCellSMDensity(kFALSE),
115  fExcludeRejectedCells(kFALSE),
116  fPlotFineGrainedCentrality(kFALSE),
117  fMaxPt(200),
118  fNCentHistBins(0),
119  fCentHistBins(0),
120  fNPtHistBins(0),
121  fPtHistBins(0),
122  fNM02HistBins(0),
123  fM02HistBins(0),
124  fUseAliEventCuts(kTRUE),
125  fEventCuts(0),
126  fEventCutList(0),
127  fUseManualEventCuts(kFALSE),
128  fGeneratorLevel(0),
129  fPHOSGeo(nullptr),
130  fHistManager(name)
131 {
133 }
134 
139 {
140 }
141 
146 {
147  fNCentHistBins = 4;
149  fCentHistBins[0] = 0;
150  fCentHistBins[1] = 10;
151  fCentHistBins[2] = 30;
152  fCentHistBins[3] = 50;
153  fCentHistBins[4] = 90;
154 
155  fNPtHistBins = 82;
158  GenerateFixedBinArray(7, 0.3, 1, fPtHistBins+6);
159  GenerateFixedBinArray(10, 1, 3, fPtHistBins+13);
160  GenerateFixedBinArray(14, 3, 10, fPtHistBins+23);
161  GenerateFixedBinArray(10, 10, 20, fPtHistBins+37);
162  GenerateFixedBinArray(15, 20, 50, fPtHistBins+47);
163  GenerateFixedBinArray(20, 50, 150, fPtHistBins+62);
164 
165  fNM02HistBins = 81;
167  GenerateFixedBinArray(35, 0, 0.7, fM02HistBins);
168  GenerateFixedBinArray(6, 0.7, 1., fM02HistBins+35);
169  GenerateFixedBinArray(20, 1., 3., fM02HistBins+41);
170  GenerateFixedBinArray(10, 3., 5., fM02HistBins+61);
171  GenerateFixedBinArray(10, 5., 10., fM02HistBins+71);
172 }
173 
179 {
181 
182  fGeneratorLevel = GetMCParticleContainer("mcparticles");
183 
185 
186  TIter next(fHistManager.GetListOfHistograms());
187  TObject* obj = 0;
188  while ((obj = next())) {
189  fOutput->Add(obj);
190  }
191 
192  // Intialize AliEventCuts
193  if (fUseAliEventCuts) {
194  fEventCutList = new TList();
195  fEventCutList ->SetOwner();
196  fEventCutList ->SetName("EventCutOutput");
197 
198  fEventCuts.OverrideAutomaticTriggerSelection(fOffTrigger);
199  if(fUseManualEventCuts==1)
200  {
201  fEventCuts.SetManualMode();
202  // Configure manual settings here
203  // ...
204  }
205  fEventCuts.AddQAplotsToList(fEventCutList);
206  fOutput->Add(fEventCutList);
207  }
208 
209  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
210 }
211 
212 /*
213  * This function allocates the histograms for the calorimeter performance study.
214  */
216 {
217 
220  }
221 
222  if (fPlotCellHistograms) {
224  }
225 
226  if (fPlotNeutralJets) {
228  }
229 
230  if (fPlotClustersInJets) {
232  }
233 
234 }
235 
236 /*
237  * This function allocates the histograms for the calorimeter performance study.
238  */
240 {
241  TString histname;
242  TString htitle;
243 
244  Double_t* clusType = new Double_t[3+1];
245  GenerateFixedBinArray(3, -0.5, 2.5, clusType);
246  const Int_t nRejBins = 32;
247  Double_t* rejReasonBins = new Double_t[nRejBins+1];
248  GenerateFixedBinArray(nRejBins, 0, nRejBins, rejReasonBins);
249  const Int_t nExBins = 200;
250  Double_t* exBins = new Double_t[nExBins+1];
251  GenerateFixedBinArray(nExBins, 0, 1, exBins);
252 
253  AliEmcalContainer* cont = 0;
254  TIter nextClusColl(&fClusterCollArray);
255  while ((cont = static_cast<AliEmcalContainer*>(nextClusColl()))) {
256 
257  // rejection reason plot, to make efficiency correction
258  histname = TString::Format("%s/hClusterRejectionReasonEMCal", cont->GetArrayName().Data());
259  htitle = histname + ";Rejection reason;#it{E}_{clus} (GeV/)";
260  TH2* hist = fHistManager.CreateTH2(histname.Data(), htitle.Data(), nRejBins, rejReasonBins, fNPtHistBins, fPtHistBins);
261  SetRejectionReasonLabels(hist->GetXaxis());
262 
263  histname = TString::Format("%s/hClusterRejectionReasonPHOS", cont->GetArrayName().Data());
264  htitle = histname + ";Rejection reason;#it{E}_{clus} (GeV/)";
265  TH2* histPhos = fHistManager.CreateTH2(histname.Data(), htitle.Data(), nRejBins, rejReasonBins, fNPtHistBins, fPtHistBins);
266  SetRejectionReasonLabels(histPhos->GetXaxis());
267 
268  histname = TString::Format("%s/hClusterRejectionReasonMC", cont->GetArrayName().Data());
269  htitle = histname + ";Rejection reason;#it{E}_{clus} (GeV/)";
270  TH2* histMC = fHistManager.CreateTH2(histname.Data(), htitle.Data(), nRejBins, rejReasonBins, fNPtHistBins, fPtHistBins);
271  SetRejectionReasonLabels(histMC->GetXaxis());
272 
273  // plot by SM
274  const Int_t nEmcalSM = 20;
275  for (Int_t sm = 0; sm < nEmcalSM; sm++) {
276  histname = TString::Format("%s/BySM/hEmcalClusEnergy_SM%d", cont->GetArrayName().Data(), sm);
277  htitle = histname + ";#it{E}_{cluster} (GeV);counts";
278  fHistManager.CreateTH1(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins);
279  }
280 
281  for (Int_t sm = 1; sm < 5; sm++) {
282  histname = TString::Format("%s/BySM/hPhosClusEnergy_SM%d", cont->GetArrayName().Data(), sm);
283  htitle = histname + ";#it{E}_{cluster} (GeV);counts";
284  fHistManager.CreateTH1(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins);
285  }
286 
287  // Plot cluster THnSparse (centrality, eta, phi, E, E-hadcorr, has matched track, M02, Ncells, passed dispersion cut)
288  Int_t dim = 0;
289  TString title[20];
290  Int_t nbins[20] = {0};
291  Double_t min[30] = {0.};
292  Double_t max[30] = {0.};
293  Double_t *binEdges[20] = {0};
294 
296  title[dim] = "Centrality %";
298  nbins[dim] = 18;
299  min[dim] = 0;
300  max[dim] = 90;
301  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
302  }
303  else {
304  nbins[dim] = fNCentHistBins;
305  binEdges[dim] = fCentHistBins;
306  min[dim] = fCentHistBins[0];
307  max[dim] = fCentHistBins[fNCentHistBins];
308  }
309  dim++;
310  }
311 
312  title[dim] = "#eta";
313  nbins[dim] = 28;
315  nbins[dim] = 100;
316  }
317  min[dim] = -0.7;
318  max[dim] = 0.7;
319  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
320  dim++;
321 
322  title[dim] = "#phi";
323  nbins[dim] = 100;
325  nbins[dim] = 357;
326  }
327  min[dim] = 1.;
328  max[dim] = 6.;
329  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
330  dim++;
331 
332  title[dim] = "#it{E}_{clus} (GeV)";
333  nbins[dim] = fNPtHistBins;
334  binEdges[dim] = fPtHistBins;
335  min[dim] = fPtHistBins[0];
336  max[dim] = fPtHistBins[fNPtHistBins];
337  dim++;
338 
340  title[dim] = "#it{E}_{clus, hadcorr} or #it{E}_{core} (GeV)";
341  nbins[dim] = fNPtHistBins;
342  binEdges[dim] = fPtHistBins;
343  min[dim] = fPtHistBins[0];
344  max[dim] = fPtHistBins[fNPtHistBins];
345  dim++;
346 
347  title[dim] = "Matched track";
348  nbins[dim] = 2;
349  min[dim] = -0.5;
350  max[dim] = 1.5;
351  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
352  dim++;
353 
354  title[dim] = "M02";
355  nbins[dim] = fNM02HistBins;
356  binEdges[dim] = fM02HistBins;
357  min[dim] = fM02HistBins[0];
358  max[dim] = fM02HistBins[fNM02HistBins];
359  dim++;
360 
361  title[dim] = "Ncells";
362  nbins[dim] = 30;
363  min[dim] = -0.5;
364  max[dim] = 29.5;
365  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
366  dim++;
367 
368  title[dim] = "Dispersion cut";
369  nbins[dim] = 2;
370  min[dim] = -0.5;
371  max[dim] = 1.5;
372  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
373  dim++;
374 
375  if (fGeneratorLevel) {
376  title[dim] = "Particle type1";
377  nbins[dim] = 8;
378  min[dim] = -0.5;
379  max[dim] = 7.5;
380  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
381  dim++;
382 
383  title[dim] = "Particle type2";
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  }
391 
392  if (fPlotEvenOddEta) {
393  title[dim] = "Even/odd eta";
394  nbins[dim] = 2;
395  min[dim] = -0.5;
396  max[dim] = 1.5;
397  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
398  dim++;
399  }
400 
402  title[dim] = "#DeltaR_{NN}";
403  nbins[dim] = 100;
404  min[dim] = 0.;
405  max[dim] = 1.;
406  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
407  dim++;
408  }
409 
410  if (fPlotClusterCone) {
411 
412  title[dim] = "Cone type";
413  nbins[dim] = 2;
414  min[dim] = -0.5;
415  max[dim] = 1.5;
416  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
417  dim++;
418 
419  title[dim] = "R";
420  nbins[dim] = 2;
421  min[dim] = 0;
422  max[dim] = 0.15;
423  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
424  dim++;
425 
426  title[dim] = "#it{E}_{cone} (GeV)";
427  nbins[dim] = fNPtHistBins;
428  binEdges[dim] = fPtHistBins;
429  min[dim] = fPtHistBins[0];
430  max[dim] = fPtHistBins[fNPtHistBins];
431  dim++;
432 
433  }
434 
435  if (fPlotCaloCentrality) {
436 
437  title[dim] = "#it{#rho}_{cell cone} (GeV)";
438  nbins[dim] = 100;
439  min[dim] = 0.;
440  max[dim] = 1000.;
441  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
442  dim++;
443 
444  title[dim] = "Ncells cone";
445  nbins[dim] = 100;
446  min[dim] = -0.5;
447  max[dim] = 99.5;
448  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
449  dim++;
450 
451  if (fPlotCellSMDensity) {
452  title[dim] = "#it{E}_{cell SM} (GeV)";
453  nbins[dim] = fNPtHistBins;
454  binEdges[dim] = fPtHistBins;
455  min[dim] = fPtHistBins[0];
456  max[dim] = fPtHistBins[fNPtHistBins];
457  dim++;
458 
459  title[dim] = "Ncells SM";
460  nbins[dim] = 100;
461  min[dim] = -0.5;
462  max[dim] = 999.5;
463  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
464  dim++;
465  }
466 
467  }
468 
469  TString thnname = TString::Format("%s/clusterObservables", cont->GetArrayName().Data());
470  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
471  for (Int_t i = 0; i < dim; i++) {
472  hn->GetAxis(i)->SetTitle(title[i]);
473  hn->SetBinEdges(i, binEdges[i]);
474  }
475 
476  // Plot Fcross distribution
477  if (fPlotExotics) {
478  histname = TString::Format("%s/hFcrossEMCal", cont->GetArrayName().Data());
479  htitle = histname + ";Centrality (%);Fcross;#it{E}_{clus} (GeV/)";
480  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, nExBins, exBins, fNPtHistBins, fPtHistBins);
481  }
482 
483  }
484 
485 }
486 
487 /*
488  * This function allocates the histograms for calo cells.
489  */
491 {
492  TString histname;
493  TString htitle;
494 
495  Double_t* clusType = new Double_t[3+1];
496  GenerateFixedBinArray(3, -0.5, 2.5, clusType);
497 
498  // centrality vs. cell energy vs. cell type (for all cells)
499  histname = TString::Format("Cells/hCellEnergyAll");
500  htitle = histname + ";#it{E}_{cell} (GeV);Centrality (%); Cluster type";
501  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins, 3, clusType);
502 
503  // centrality vs. cell energy vs. cell type (for cells in accepted clusters)
504  histname = TString::Format("Cells/hCellEnergyAccepted");
505  htitle = histname + ";#it{E}_{cell} (GeV);Centrality (%); Cluster type";
506  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins, 3, clusType);
507 
508  // centrality vs. cell energy vs. cell type (for leading cells in accepted clusters)
509  histname = TString::Format("Cells/hCellEnergyLeading");
510  htitle = histname + ";#it{E}_{cell} (GeV);Centrality (%); Cluster type";
511  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins, 3, clusType);
512 
513  // plot cell patches by SM
514  const Int_t nEmcalSM = 20;
515  for (Int_t sm = 0; sm < nEmcalSM; sm++) {
516  histname = TString::Format("Cells/BySM/hEmcalPatchEnergy_SM%d", sm);
517  htitle = histname + ";#it{E}_{cell patch} (GeV);Centrality (%)";
518  fHistManager.CreateTH2(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins);
519  }
520 
521  for (Int_t sm = 1; sm < 5; sm++) {
522  histname = TString::Format("Cells/BySM/hPhosPatchEnergy_SM%d", sm);
523  htitle = histname + ";#it{E}_{cell patch} (GeV);Centrality (%)";
524  fHistManager.CreateTH2(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins);
525  }
526 
527 }
528 
529 /*
530  * This function allocates the histograms for neutral jets.
531  */
533 {
534  AliJetContainer* jets = 0;
535  TIter nextJetColl(&fJetCollArray);
536  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
537 
538  TString axisTitle[30]= {""};
539  Int_t nbinsJet[30] = {0};
540  Double_t minJet[30] = {0.};
541  Double_t maxJet[30] = {0.};
542  Double_t *binEdgesJet[20] = {0};
543  Int_t dimJet = 0;
544 
545  if (fForceBeamType != kpp) {
546  axisTitle[dimJet] = "Centrality (%)";
547  nbinsJet[dimJet] = fNCentHistBins;
548  binEdgesJet[dimJet] = fCentHistBins;
549  minJet[dimJet] = fCentHistBins[0];
550  maxJet[dimJet] = fCentHistBins[fNCentHistBins];
551  dimJet++;
552  }
553 
554  axisTitle[dimJet] = "#eta_{jet}";
555  nbinsJet[dimJet] = 28;
556  minJet[dimJet] = -0.7;
557  maxJet[dimJet] = 0.7;
558  binEdgesJet[dimJet] = GenerateFixedBinArray(nbinsJet[dimJet], minJet[dimJet], maxJet[dimJet]);
559  dimJet++;
560 
561  axisTitle[dimJet] = "#phi_{jet} (rad)";
562  nbinsJet[dimJet] = 100;
563  minJet[dimJet] = 1.;
564  maxJet[dimJet] = 6.;
565  binEdgesJet[dimJet] = GenerateFixedBinArray(nbinsJet[dimJet], minJet[dimJet], maxJet[dimJet]);
566  dimJet++;
567 
568  axisTitle[dimJet] = "#it{E}_{T} (GeV)";
569  nbinsJet[dimJet] = fNPtHistBins;
570  binEdgesJet[dimJet] = fPtHistBins;
571  minJet[dimJet] = fPtHistBins[0];
572  maxJet[dimJet] = fPtHistBins[fNPtHistBins];
573  dimJet++;
574 
575  axisTitle[dimJet] = "#rho (GeV/#it{c})";
576  nbinsJet[dimJet] = 100;
577  minJet[dimJet] = 0.;
578  maxJet[dimJet] = 1000.;
579  binEdgesJet[dimJet] = GenerateFixedBinArray(nbinsJet[dimJet], minJet[dimJet], maxJet[dimJet]);
580  dimJet++;
581 
582  axisTitle[dimJet] = "N_{clusters}";
583  nbinsJet[dimJet] = 20;
584  minJet[dimJet] = -0.5;
585  maxJet[dimJet] = 19.5;
586  binEdgesJet[dimJet] = GenerateFixedBinArray(nbinsJet[dimJet], minJet[dimJet], maxJet[dimJet]);
587  dimJet++;
588 
589  axisTitle[dimJet] = "#it{E}_{T}, acc clus within R (GeV)";
590  nbinsJet[dimJet] = fNPtHistBins;
591  binEdgesJet[dimJet] = fPtHistBins;
592  minJet[dimJet] = fPtHistBins[0];
593  maxJet[dimJet] = fPtHistBins[fNPtHistBins];
594  dimJet++;
595 
596  axisTitle[dimJet] = "#it{E}_{T}, acc cell 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  TString thnname = TString::Format("%s/hClusterJetObservables", jets->GetArrayName().Data());
604  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dimJet, nbinsJet, minJet, maxJet);
605  for (Int_t i = 0; i < dimJet; i++) {
606  hn->GetAxis(i)->SetTitle(axisTitle[i]);
607  hn->SetBinEdges(i, binEdgesJet[i]);
608  }
609 
610  }
611 }
612 
613 /*
614  * This function allocates the histograms for clusters within jets.
615  */
617 {
618  TString histname;
619  TString htitle;
620  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
621 
622  AliJetContainer* jets = 0;
623  TIter nextJetColl(&fJetCollArray);
624  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
625 
626  // Plot cluster spectra within jets
627  // (centrality, cluster energy, jet pT, jet eta, jet phi)
628  Int_t dim = 0;
629  TString title[20];
630  Int_t nbins[20] = {0};
631  Double_t min[30] = {0.};
632  Double_t max[30] = {0.};
633  Double_t *binEdges[20] = {0};
634 
636  title[dim] = "Centrality %";
637  nbins[dim] = fNCentHistBins;
638  binEdges[dim] = fCentHistBins;
639  min[dim] = fCentHistBins[0];
640  max[dim] = fCentHistBins[fNCentHistBins];
641  dim++;
642  }
643 
644  title[dim] = "#it{E}_{clus} (GeV)";
645  nbins[dim] = fNPtHistBins;
646  binEdges[dim] = fPtHistBins;
647  min[dim] = fPtHistBins[0];
648  max[dim] = fPtHistBins[fNPtHistBins];
649  dim++;
650 
651  title[dim] = "#it{p}_{T,jet}^{corr}";
652  nbins[dim] = nPtBins;
653  min[dim] = 0;
654  max[dim] = fMaxPt;
655  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
656  dim++;
657 
658  title[dim] = "#eta_{jet}";
659  nbins[dim] = 40;
660  min[dim] = -0.5;
661  max[dim] = 0.5;
662  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
663  dim++;
664 
665  title[dim] = "#phi_{jet}";
666  nbins[dim] = 200;
667  min[dim] = 1.;
668  max[dim] = 6.;
669  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
670  dim++;
671 
672  TString thnname = TString::Format("%s/hClustersInJets", jets->GetArrayName().Data());
673  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
674  for (Int_t i = 0; i < dim; i++) {
675  hn->GetAxis(i)->SetTitle(title[i]);
676  hn->SetBinEdges(i, binEdges[i]);
677  }
678 
679  // (jet type, jet pT, cluster shift)
680  histname = TString::Format("%s/hCaloJESshift", jets->GetArrayName().Data());
681  htitle = histname + ";type;#it{p}_{T}^{corr} (GeV/#it{c});#Delta_{JES}";
682  fHistManager.CreateTH3(histname.Data(), htitle.Data(), 3, -0.5, 2.5, nPtBins, 0, fMaxPt, 100, 0, 20);
683 
684  }
685 }
686 
692 {
694 
695  fNeedEmcalGeom = kTRUE;
696 
697  // Load the PHOS geometry
698  fPHOSGeo = AliPHOSGeometry::GetInstance();
699  if (fPHOSGeo) {
700  AliInfo("Found instance of PHOS geometry!");
701  }
702  else {
703  AliInfo("Creating PHOS geometry!");
704  Int_t runNum = InputEvent()->GetRunNumber();
705  if(runNum<209122) //Run1
706  fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP");
707  else
708  fPHOSGeo = AliPHOSGeometry::GetInstance("Run2");
709 
710  if (fPHOSGeo) {
711  AliOADBContainer geomContainer("phosGeo");
712  geomContainer.InitFromFile("$ALICE_PHYSICS/OADB/PHOS/PHOSMCGeometry.root","PHOSMCRotationMatrixes");
713  TObjArray* matrixes = (TObjArray*)geomContainer.GetObject(runNum,"PHOSRotationMatrixes");
714  for(Int_t mod=0; mod<6; mod++) {
715  if(!matrixes->At(mod)) continue;
716  fPHOSGeo->SetMisalMatrix(((TGeoHMatrix*)matrixes->At(mod)),mod);
717  printf(".........Adding Matrix(%d), geo=%p\n",mod,fPHOSGeo);
718  ((TGeoHMatrix*)matrixes->At(mod))->Print();
719  }
720  }
721  }
722 }
723 
728 {
729  if (fUseAliEventCuts) {
730  if (!fEventCuts.AcceptEvent(InputEvent()))
731  {
732  PostData(1, fOutput);
733  return kFALSE;
734  }
735  }
736  else {
738  }
739  return kTRUE;
740 }
741 
750 {
751  return kTRUE;
752 }
753 
761 {
762 
765  }
766 
767  if (fPlotNeutralJets) {
769  }
770 
771  if (fPlotClustersInJets) {
773  }
774 
775  if (fPlotCellHistograms) {
777  }
778 
779  return kTRUE;
780 }
781 
782 /*
783  * This function fills the histograms for the calorimeter performance study.
784  */
786 {
787  // Define some vars
788  TString histname;
789  Double_t Enonlin;
790  Double_t Ehadcorr;
791  Int_t absId;
792  Double_t ecell;
793  Double_t leadEcell;
794  Int_t leadAbsId;
795 
796  // Get cells from event
797  fCaloCells = InputEvent()->GetEMCALCells();
798  AliVCaloCells* phosCaloCells = InputEvent()->GetPHOSCells();
799 
800  // If MC, get the MC event
801  const AliMCEvent* mcevent = nullptr;
802  if (fGeneratorLevel) {
803  mcevent = MCEvent();
804  }
805 
806  // Loop through clusters and plot cluster THnSparse (centrality, cluster type, E, E-hadcorr, has matched track, M02, Ncells)
807  AliClusterContainer* clusters = 0;
808  const AliVCluster* clus;
809  TString clustersName;
810  TIter nextClusColl(&fClusterCollArray);
811  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
813  clustersName = clusters->GetArrayName();
814  for (AliClusterIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
815 
816  clus = it->second;
817 
818  // Determine cluster type (EMCal/DCal/Phos)
819  ClusterType clusType = kNA;
820  if (clus->IsEMCAL()) {
821  Double_t phi = it->first.Phi_0_2pi();
822  Int_t isDcal = Int_t(phi > fgkEMCalDCalPhiDivide);
823  if (isDcal == 0) {
824  clusType = kEMCal;
825  } else if (isDcal == 1) {
826  clusType = kDCal;
827  }
828  } else if (clus->GetType() == AliVCluster::kPHOSNeutral){
829  clusType = kPHOS;
830  }
831 
832  // rejection reason plots, to make efficiency correction
833  if (clus->IsEMCAL()) {
834  histname = TString::Format("%s/hClusterRejectionReasonEMCal", clusters->GetArrayName().Data());
835  UInt_t rejectionReason = 0;
836  if (!clusters->AcceptCluster(it.current_index(), rejectionReason)) {
837  fHistManager.FillTH2(histname, clusters->GetRejectionReasonBitPosition(rejectionReason), it->first.E());
838  continue;
839  }
840  } else if (clus->GetType() == AliVCluster::kPHOSNeutral){
841  histname = TString::Format("%s/hClusterRejectionReasonPHOS", clusters->GetArrayName().Data());
842  UInt_t rejectionReason = 0;
843  if (!clusters->AcceptCluster(it.current_index(), rejectionReason)) {
844  fHistManager.FillTH2(histname, clusters->GetRejectionReasonBitPosition(rejectionReason), it->first.E());
845  continue;
846  }
847  } else {
848  continue; // avoid CPV clusters
849  }
850 
851  // Fill cluster spectra by SM, and fill cell histograms
852  Enonlin = 0;
853  Ehadcorr = 0;
854  if (clus->IsEMCAL()) {
855 
856  Ehadcorr = clus->GetHadCorrEnergy();
857  Enonlin = clus->GetNonLinCorrEnergy();
859  Enonlin = clus->E();
860  }
861 
862  if (fPlotExotics) {
863  histname = TString::Format("%s/hFcrossEMCal", clusters->GetArrayName().Data());
864  Double_t Fcross = GetFcross(clus, fCaloCells);
865  fHistManager.FillTH3(histname, fCent, Fcross, Enonlin);
866  }
867 
868  Int_t sm = fGeom->GetSuperModuleNumber(clus->GetCellAbsId(0));
869  if (sm >=0 && sm < 20) {
870  histname = TString::Format("%s/BySM/hEmcalClusEnergy_SM%d", clusters->GetArrayName().Data(), sm);
871  fHistManager.FillTH1(histname, Enonlin);
872  }
873  else {
874  AliError(Form("Supermodule %d does not exist!", sm));
875  }
876 
877  // Get cells from each accepted cluster, and plot centrality vs. cell energy vs. cell type
878  histname = TString::Format("Cells/hCellEnergyAccepted");
879  leadEcell = 0;
880  for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++){
881  absId = clus->GetCellAbsId(iCell);
882  ecell = fCaloCells->GetCellAmplitude(absId);
883  fHistManager.FillTH3(histname, ecell, fCent, kEMCal); // Note: I don't distinguish EMCal from DCal cells
884  if (ecell > leadEcell) {
885  leadEcell = ecell;
886  leadAbsId = absId;
887  }
888  }
889  // Plot also the leading cell
890  histname = TString::Format("Cells/hCellEnergyLeading");
891  fHistManager.FillTH3(histname, leadEcell, fCent, kEMCal);
892 
893  } else if (clus->GetType() == AliVCluster::kPHOSNeutral){
894 
895  Ehadcorr = clus->GetCoreEnergy();
896  Enonlin = clus->E();
897 
898  Int_t relid[4];
899  if (fPHOSGeo) {
900  fPHOSGeo->AbsToRelNumbering(clus->GetCellAbsId(0), relid);
901  Int_t sm = relid[0];
902  if (sm >=1 && sm < 5) {
903  histname = TString::Format("%s/BySM/hPhosClusEnergy_SM%d", clusters->GetArrayName().Data(), sm);
904  fHistManager.FillTH1(histname, clus->E());
905  }
906  else {
907  AliError(Form("Supermodule %d does not exist!", sm));
908  }
909  }
910 
911  // Get cells from each accepted cluster, and plot centrality vs. cell energy vs. cell type
912  histname = TString::Format("Cells/hCellEnergyAccepted");
913  leadEcell = 0;
914  for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++){
915  absId = clus->GetCellAbsId(iCell);
916  ecell = phosCaloCells->GetCellAmplitude(absId);
917  fHistManager.FillTH3(histname, ecell, fCent, kPHOS);
918  if (ecell > leadEcell) {
919  leadEcell = ecell;
920  }
921  }
922  // Plot also the leading cell
923  histname = TString::Format("Cells/hCellEnergyLeading");
924  fHistManager.FillTH3(histname, leadEcell, fCent, kPHOS);
925  }
926 
927  // Check if the cluster has a matched track
928  Int_t hasMatchedTrack = -1;
929  Int_t nMatchedTracks = clus->GetNTracksMatched();
930  if (nMatchedTracks == 0) {
931  hasMatchedTrack = 0;
932  } else if (nMatchedTracks > 0) {
933  hasMatchedTrack = 1;
934  }
935 
936  // Check if the cluster passes the dispersion cut for photon-like cluster (meaningful only for PHOS)
937  Int_t passedDispersionCut = 0;
938  if (clus->Chi2() < 2.5*2.5) {
939  passedDispersionCut = 1;
940  }
941 
942  // Fill info about the cluster
943  Double_t eta = it->first.Eta();
944  Double_t phi = it->first.Phi_0_2pi();
945  Double_t M02 = clus->GetM02();
946  Int_t nCells = clus->GetNCells();
947  Double_t distNN = FindNearestNeighborDistance(it->first);
948 
949  // If cluster is EMCal, find whether the eta column is even or odd
950  Int_t isOddEta = -1;
951  Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta;
952  if (clus->IsEMCAL()) {
953  fGeom->GetCellIndex(leadAbsId, nSupMod, nModule, nIphi, nIeta);
954  fGeom->GetCellPhiEtaIndexInSModule(nSupMod, nModule, nIphi, nIeta, iphi, ieta);
955  isOddEta = ieta % 2;
956  }
957 
958  // Standard option: fill once per cluster
960  FillClusterTHnSparse(clustersName, eta, phi, Enonlin, Ehadcorr, hasMatchedTrack, M02, nCells, passedDispersionCut, distNN, isOddEta);
961  }
962 
963  if (fPlotCaloCentrality) {
964 
965  Double_t eCellCone = 0.;
966  Int_t nCellsCone = 0;
967  Double_t eCellSM = 0.;
968  Int_t nCellsSM = 0;
969  Double_t areaCone = TMath::Pi() * 0.07 * 0.07;
970  Double_t areaCell;
971  Double_t eDensityCone;
972 
973  // Get the SM number
974  Int_t sm = -1;
975  if (clusType == kEMCal) {
976  sm = fGeom->GetSuperModuleNumber(clus->GetCellAbsId(0));
977  areaCell = 0.014*0.014;
978  }
979  if (clusType == kPHOS) {
980  Int_t relid[4];
981  fPHOSGeo->AbsToRelNumbering(clus->GetCellAbsId(0), relid);
982  sm = relid[0];
983  areaCell = 0.014*0.014*(2.2/6.0)*(2.2/6.0); // approximating same r
984  }
985 
986  // Only fill the THnSparse if the cluster is located in a full SM of EMCal or PHOS
987  if ( (clusType == kEMCal && sm < 10 ) || (clusType == kPHOS && sm < 4) ) {
988 
989  eCellCone = GetConeCellEnergy(eta, phi, 0.07) - Enonlin;
990  nCellsCone = (Int_t)GetConeCellEnergy(eta, phi, 0.07, kTRUE) - nCells;
991 
992  eDensityCone = eCellCone / (areaCone - nCells*areaCell);
993 
994  if (fPlotCellSMDensity) {
995  eCellSM = GetSMCellEnergy(sm, clusType) - Enonlin;
996  nCellsSM = (Int_t)GetSMCellEnergy(sm, clusType, kTRUE) - nCells;
997  }
998 
999  FillClusterTHnSparse(clustersName, eta, phi, Enonlin, eDensityCone, eCellSM, nCellsCone, nCellsSM);
1000 
1001  }
1002 
1003  }
1004 
1005  // If cluster cone option enabled, fill for each R and cone type
1006  if (fPlotClusterCone) {
1007 
1008  // cluster cone, R=0.05
1009  FillClusterTHnSparse(clustersName, eta, phi, Enonlin, Ehadcorr, hasMatchedTrack, M02, nCells, passedDispersionCut, distNN, isOddEta, 0, 0.05, GetConeClusterEnergy(eta, phi, 0.05));
1010  // cluster cone, R=0.1
1011  FillClusterTHnSparse(clustersName, eta, phi, Enonlin, Ehadcorr, hasMatchedTrack, M02, nCells, passedDispersionCut, distNN, isOddEta, 0, 0.1, GetConeClusterEnergy(eta, phi, 0.1));
1012  // cell cone, R=0.05
1013  FillClusterTHnSparse(clustersName, eta, phi, Enonlin, Ehadcorr, hasMatchedTrack, M02, nCells, passedDispersionCut, distNN, isOddEta, 1, 0.05, GetConeCellEnergy(eta, phi, 0.05));
1014  // cell cone, R=0.1
1015  FillClusterTHnSparse(clustersName, eta, phi, Enonlin, Ehadcorr, hasMatchedTrack, M02, nCells, passedDispersionCut, distNN, isOddEta, 1, 0.1, GetConeCellEnergy(eta, phi, 0.1));
1016 
1017  }
1018 
1019  }
1020  }
1021 }
1022 
1023 /*
1024  * This function fills the cluster THnSparse.
1025  */
1026 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 coneType, Double_t R, Double_t Econe)
1027 {
1028  Double_t contents[30]={0};
1029  TString histname = TString::Format("%s/clusterObservables", clustersName.Data());
1030  THnSparse* histClusterObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1031  if (!histClusterObservables) return;
1032  for (Int_t i = 0; i < histClusterObservables->GetNdimensions(); i++) {
1033  TString title(histClusterObservables->GetAxis(i)->GetTitle());
1034  if (title=="Centrality %")
1035  contents[i] = fCent;
1036  else if (title=="#eta")
1037  contents[i] = eta;
1038  else if (title=="#phi")
1039  contents[i] = phi;
1040  else if (title=="#it{E}_{clus} (GeV)")
1041  contents[i] = Enonlin;
1042  else if (title=="#it{E}_{clus, hadcorr} or #it{E}_{core} (GeV)")
1043  contents[i] = Ehadcorr;
1044  else if (title=="Matched track")
1045  contents[i] = hasMatchedTrack;
1046  else if (title=="M02")
1047  contents[i] = M02;
1048  else if (title=="Ncells")
1049  contents[i] = nCells;
1050  else if (title=="Dispersion cut")
1051  contents[i] = passedDispersionCut;
1052  else if (title=="#DeltaR_{NN}")
1053  contents[i] = distNN;
1054  else if (title=="Even/odd eta")
1055  contents[i] = isOddEta;
1056  else if (title=="Cone type")
1057  contents[i] = coneType;
1058  else if (title=="R")
1059  contents[i] = R;
1060  else if (title=="#it{E}_{cone} (GeV)")
1061  contents[i] = Econe;
1062  else
1063  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1064  }
1065  histClusterObservables->Fill(contents);
1066 
1067 }
1068 
1069 /*
1070  * This function fills the cluster THnSparse (alternate signature, used for local density option).
1071  */
1072 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)
1073 {
1074  Double_t contents[30]={0};
1075  TString histname = TString::Format("%s/clusterObservables", clustersName.Data());
1076  THnSparse* histClusterObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1077  if (!histClusterObservables) return;
1078  for (Int_t i = 0; i < histClusterObservables->GetNdimensions(); i++) {
1079  TString title(histClusterObservables->GetAxis(i)->GetTitle());
1080  if (title=="Centrality %")
1081  contents[i] = fCent;
1082  else if (title=="#eta")
1083  contents[i] = eta;
1084  else if (title=="#phi")
1085  contents[i] = phi;
1086  else if (title=="#it{E}_{clus} (GeV)")
1087  contents[i] = Enonlin;
1088  else if (title=="#it{#rho}_{cell cone} (GeV)")
1089  contents[i] = eCellCone;
1090  else if (title=="Ncells cone")
1091  contents[i] = nCellsCone;
1092  else if (title=="#it{E}_{cell SM} (GeV)")
1093  contents[i] = eCellSM;
1094  else if (title=="Ncells SM")
1095  contents[i] = nCellsSM;
1096  else
1097  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1098  }
1099  histClusterObservables->Fill(contents);
1100 
1101 }
1102 
1103 /*
1104  * This function fills cell histograms.
1105  */
1107 {
1108  // Define some vars
1109  TString histname;
1110  Int_t absId;
1111  Double_t ecell;
1112 
1113  // Get cells from event
1114  fCaloCells = InputEvent()->GetEMCALCells();
1115  AliVCaloCells* phosCaloCells = InputEvent()->GetPHOSCells();
1116 
1117  // Loop through all cells and fill histos
1118  Int_t sm;
1119  Int_t relid[4];
1120  Double_t patchSumEMCal[20] = {0.};
1121  Double_t patchSumPHOS[4] = {0.};
1122  for (Int_t i=0; i<fCaloCells->GetNumberOfCells(); i++) {
1123 
1124  absId = fCaloCells->GetCellNumber(i);
1125  ecell = fCaloCells->GetCellAmplitude(absId);
1126 
1127  // Fill cell histo
1128  histname = TString::Format("Cells/hCellEnergyAll");
1129  fHistManager.FillTH3(histname, ecell, fCent, kEMCal); // Note: I don't distinguish EMCal from DCal cells
1130 
1131  // Fill cell patch histo, per SM
1132  sm = fGeom->GetSuperModuleNumber(absId);
1133  if (sm >=0 && sm < 20) {
1134  patchSumEMCal[sm] += ecell;
1135  }
1136 
1137  }
1138 
1139  for (Int_t i=0; i<phosCaloCells->GetNumberOfCells(); i++) {
1140 
1141  absId = phosCaloCells->GetCellNumber(i);
1142  ecell = phosCaloCells->GetCellAmplitude(absId);
1143 
1144  if (absId < 0) {
1145  continue; // skip CPV cells
1146  }
1147 
1148  // Fill cell histo
1149  histname = TString::Format("Cells/hCellEnergyAll");
1150  fHistManager.FillTH3(histname, ecell, fCent, kPHOS);
1151 
1152  // Fill cell patch histo, per SM
1153  fPHOSGeo->AbsToRelNumbering(absId, relid);
1154  sm = relid[0];
1155  if (sm >=1 && sm < 5) {
1156  patchSumPHOS[sm-1] += ecell;
1157  }
1158 
1159  }
1160 
1161  for (Int_t sm = 0; sm < 20; sm++) {
1162  histname = TString::Format("Cells/BySM/hEmcalPatchEnergy_SM%d", sm);
1163  fHistManager.FillTH2(histname, patchSumEMCal[sm], fCent);
1164  }
1165 
1166  for (Int_t sm = 1; sm < 5; sm++) {
1167  histname = TString::Format("Cells/BySM/hPhosPatchEnergy_SM%d", sm);
1168  fHistManager.FillTH2(histname, patchSumPHOS[sm-1], fCent);
1169  }
1170 
1171 }
1172 
1173 /*
1174  * This function fills neutral jet histograms.
1175  */
1177 {
1178  AliJetContainer* jets = 0;
1179  TIter nextJetColl(&fJetCollArray);
1180  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
1181 
1182  // plot neutral jets THnSparse (centrality, eta, phi, E, Nclusters)
1183  Double_t contents[30]={0};
1184  TString histname = TString::Format("%s/hClusterJetObservables", jets->GetArrayName().Data());
1185  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1186  if (!histJetObservables) return;
1187 
1188  for (auto jet : jets->accepted()) {
1189 
1190  for (Int_t i = 0; i < histJetObservables->GetNdimensions(); i++) {
1191  TString title(histJetObservables->GetAxis(i)->GetTitle());
1192  if (title=="Centrality (%)")
1193  contents[i] = fCent;
1194  else if (title=="#eta_{jet}")
1195  contents[i] = jet->Eta();
1196  else if (title=="#phi_{jet} (rad)")
1197  contents[i] = jet->Phi_0_2pi();
1198  else if (title=="#it{E}_{T} (GeV)")
1199  contents[i] = jet->Pt();
1200  else if (title=="#rho (GeV/#it{c})")
1201  contents[i] = jet->Pt() / jet->Area();
1202  else if (title=="N_{clusters}")
1203  contents[i] = jet->GetNumberOfClusters();
1204  else if (title=="#it{E}_{T}, acc clus within R (GeV)")
1205  contents[i] = GetConeClusterEnergy(jet->Eta(), jet->Phi_0_2pi(), jets->GetJetRadius());
1206  else if (title=="#it{E}_{T}, acc cell within R (GeV)")
1207  contents[i] = GetConeCellEnergy(jet->Eta(), jet->Phi_0_2pi(), jets->GetJetRadius());
1208  else
1209  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1210  }
1211  histJetObservables->Fill(contents);
1212 
1213  }
1214  }
1215 }
1216 
1217 /*
1218  * This function fills clusters within jets and estimates the JES shift due to the bump.
1219  */
1221 {
1222  TString histname;
1223  Double_t rho;
1224  AliJetContainer* jets = 0;
1225  TIter nextJetColl(&fJetCollArray);
1226  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
1227 
1228  rho = jets->GetRhoVal();
1229 
1230  for (const auto jet : jets->accepted()) {
1231 
1232  // Fill cluster spectra of clusters within jets
1233  //(centrality, cluster energy, jet pT, jet eta, jet phi)
1234  histname = TString::Format("%s/hClustersInJets", jets->GetArrayName().Data());
1235  Int_t nClusters = jet->GetNumberOfClusters();
1236  AliVCluster* clus;
1237  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
1238 
1239  clus = jet->Cluster(iClus);
1240 
1242  Double_t x[5] = {fCent, clus->GetNonLinCorrEnergy(), GetJetPt(jet, rho), jet->Eta(), jet->Phi_0_2pi()};
1243  fHistManager.FillTHnSparse(histname, x);
1244  }
1245  else {
1246  Double_t x[4] = {clus->GetNonLinCorrEnergy(), GetJetPt(jet, rho), jet->Eta(), jet->Phi_0_2pi()};
1247  fHistManager.FillTHnSparse(histname, x);
1248  }
1249 
1250  }
1251 
1252  // Loop through clusters, and plot estimated shift in JES due to cluster bump
1253  // Only do for 0-10% centrality, and for EMCal/DCal
1254  Double_t eclus;
1255  Double_t shift;
1256  Double_t shiftSum = 0;
1258  if (GetJetType(jet) > -0.5 && GetJetType(jet) < 1.5) {
1259  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
1260  clus = jet->Cluster(iClus);
1261  eclus = clus->GetNonLinCorrEnergy();
1262  if (eclus > 0.5) {
1263  shift = 0.79 * TMath::Exp(-0.5 * ((eclus - 3.81) / 1.50)*((eclus - 3.81) / 1.50) );
1264  shiftSum += shift;
1265  }
1266  }
1267  histname = TString::Format("%s/hCaloJESshift", jets->GetArrayName().Data());
1268  fHistManager.FillTH3(histname, GetJetType(jet), GetJetPt(jet, rho), shiftSum);
1269  }
1270  }
1271 
1272  }
1273  }
1274 }
1275 
1280 {
1281  Double_t pT = jet->Pt() - rho * jet->Area();
1282  return pT;
1283 }
1284 
1289 {
1290  Double_t deltaPhi = TMath::Abs(part.Phi_0_2pi() - phiRef);
1291  Double_t deltaEta = TMath::Abs(part.Eta() - etaRef);
1292  Double_t deltaR = TMath::Sqrt( deltaPhi*deltaPhi + deltaEta*deltaEta );
1293  return deltaR;
1294 }
1295 
1300 {
1301  Double_t deltaPhi = TMath::Abs(phi1-phi2);
1302  Double_t deltaEta = TMath::Abs(eta1-eta2);
1303  Double_t deltaR = TMath::Sqrt( deltaPhi*deltaPhi + deltaEta*deltaEta );
1304  return deltaR;
1305 }
1306 
1311 {
1312  UInt_t jetType = jet->GetJetAcceptanceType();
1313  Double_t type = -1;
1314  if (jetType & AliEmcalJet::kEMCAL) {
1315  type = kEMCal;
1316  }
1317  else if (jetType & AliEmcalJet::kDCALonly) {
1318  type = kDCal;
1319  }
1320  else if (jetType & AliEmcalJet::kPHOS) {
1321  type = kPHOS;
1322  }
1323 
1324  return type;
1325 }
1326 
1330 Double_t AliAnalysisTaskEmcalVsPhos::GetFcross(const AliVCluster *cluster, AliVCaloCells *cells)
1331 {
1332  Int_t AbsIdseed = -1;
1333  Double_t Eseed = 0;
1334  for (Int_t i = 0; i < cluster->GetNCells(); i++) {
1335  if (cells->GetCellAmplitude(cluster->GetCellAbsId(i)) > Eseed) {
1336  Eseed = cells->GetCellAmplitude(cluster->GetCellAbsId(i));
1337  AbsIdseed = cluster->GetCellAbsId(i);
1338  }
1339  }
1340 
1341  if (Eseed < 1e-9) {
1342  return 100;
1343  }
1344 
1345  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
1346  fGeom->GetCellIndex(AbsIdseed,imod,iTower,iIphi,iIeta);
1347  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi,iIeta,iphi,ieta);
1348 
1349  //Get close cells index and energy, not in corners
1350 
1351  Int_t absID1 = -1;
1352  Int_t absID2 = -1;
1353 
1354  if (iphi < AliEMCALGeoParams::fgkEMCALRows-1) {
1355  absID1 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
1356  }
1357  if (iphi > 0) {
1358  absID2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
1359  }
1360 
1361  // In case of cell in eta = 0 border, depending on SM shift the cross cell index
1362 
1363  Int_t absID3 = -1;
1364  Int_t absID4 = -1;
1365 
1366  if (ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2)) {
1367  absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
1368  absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
1369  }
1370  else if (ieta == 0 && imod%2) {
1371  absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
1372  absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
1373  }
1374  else {
1375  if (ieta < AliEMCALGeoParams::fgkEMCALCols-1) {
1376  absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
1377  }
1378  if (ieta > 0) {
1379  absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
1380  }
1381  }
1382 
1383  Double_t ecell1 = cells->GetCellAmplitude(absID1);
1384  Double_t ecell2 = cells->GetCellAmplitude(absID2);
1385  Double_t ecell3 = cells->GetCellAmplitude(absID3);
1386  Double_t ecell4 = cells->GetCellAmplitude(absID4);
1387 
1388  Double_t Ecross = ecell1 + ecell2 + ecell3 + ecell4;
1389 
1390  Double_t Fcross = 1 - Ecross/Eseed;
1391 
1392  return Fcross;
1393 }
1394 
1399 {
1400  Double_t distNN = 10.;
1401  Double_t etaRef = clusterRef.Eta();
1402  Double_t phiRef = clusterRef.Phi_0_2pi();
1403 
1405  AliTLorentzVector clusNNcand;
1406  for (auto clusIterator : clusters->accepted_momentum() ) {
1407 
1408  clusNNcand.Clear();
1409  clusNNcand = clusIterator.first;
1410 
1411  Double_t distNNcand = GetDeltaR(clusNNcand, etaRef, phiRef);
1412 
1413  if (distNNcand < distNN && distNNcand > 0.001) {
1414  distNN = distNNcand;
1415  }
1416 
1417  }
1418 
1419  return distNN;
1420 
1421 }
1422 
1427 {
1429  AliTLorentzVector clus;
1430  Double_t energy = 0.;
1431  for (auto clusIterator : clusCont->accepted_momentum() ) {
1432 
1433  clus.Clear();
1434  clus = clusIterator.first;
1435 
1436  if (GetDeltaR(clus, etaRef, phiRef) < R) {
1437  energy += clus.E();
1438  }
1439  }
1440  return energy;
1441 }
1442 
1448 {
1449  Double_t energy = 0.;
1450  Double_t nCells = 0.;
1451 
1452  // Get cells from event
1453  fCaloCells = InputEvent()->GetEMCALCells();
1454  AliVCaloCells* phosCaloCells = InputEvent()->GetPHOSCells();
1455 
1456  Double_t eta;
1457  Double_t phi;
1458  Int_t absId;
1459  TVector3 pos;
1460  Int_t relid[4];
1461  Int_t sm;
1462 
1463  for (Int_t i=0; i<fCaloCells->GetNumberOfCells(); i++) {
1464 
1465  absId = fCaloCells->GetCellNumber(i);
1466 
1467  fGeom->EtaPhiFromIndex(absId, eta, phi);
1468  phi = TVector2::Phi_0_2pi(phi);
1469 
1470  if (GetDeltaR(eta, phi, etaRef, phiRef) < R) {
1471 
1472  if (fExcludeRejectedCells) {
1473  if (IsCellRejected(absId, kEMCal)) {
1474  continue;
1475  }
1476  }
1477 
1478  energy += fCaloCells->GetCellAmplitude(absId);
1479  nCells += 1.;
1480  }
1481 
1482  }
1483 
1484  for (Int_t i=0; i<phosCaloCells->GetNumberOfCells(); i++) {
1485 
1486  absId = phosCaloCells->GetCellNumber(i);
1487 
1488  fPHOSGeo->AbsToRelNumbering(absId, relid);
1489  sm = relid[0];
1490  if (sm < 1 || sm > 4) {
1491  continue;
1492  }
1493 
1494  fPHOSGeo->RelPosInAlice(absId, pos); // pos then contains (x,y,z) coordinate of cell
1495  eta = pos.Eta();
1496  phi = pos.Phi();
1497  phi = TVector2::Phi_0_2pi(phi);
1498 
1499  if (GetDeltaR(eta, phi, etaRef, phiRef) < R) {
1500 
1501  if (fExcludeRejectedCells) {
1502  if (IsCellRejected(absId, kPHOS)) {
1503  continue;
1504  }
1505  }
1506 
1507  energy += phosCaloCells->GetCellAmplitude(absId);
1508  nCells += 1.;
1509  }
1510 
1511  }
1512 
1513  if (returnNcells) {
1514  return nCells;
1515  }
1516 
1517  return energy;
1518 }
1519 
1525 {
1526  Double_t energy = 0.;
1527  Double_t nCells = 0.;
1528  Int_t absId;
1529  Int_t cellSM;
1530  Int_t relid[4];
1531 
1532  if (clusType == kEMCal) {
1533 
1534  fCaloCells = InputEvent()->GetEMCALCells();
1535 
1536  for (Int_t i=0; i<fCaloCells->GetNumberOfCells(); i++) {
1537 
1538  absId = fCaloCells->GetCellNumber(i);
1539  cellSM = fGeom->GetSuperModuleNumber(absId);
1540 
1541  if (cellSM == sm) {
1542 
1543  if (fExcludeRejectedCells) {
1544  if (IsCellRejected(absId, kEMCal)) {
1545  continue;
1546  }
1547  }
1548 
1549  energy += fCaloCells->GetCellAmplitude(absId);
1550  nCells += 1.;
1551 
1552  }
1553  }
1554  }
1555 
1556  if (clusType == kPHOS) {
1557 
1558  AliVCaloCells* phosCaloCells = InputEvent()->GetPHOSCells();
1559 
1560  for (Int_t i=0; i<phosCaloCells->GetNumberOfCells(); i++) {
1561 
1562  absId = phosCaloCells->GetCellNumber(i);
1563 
1564  fPHOSGeo->AbsToRelNumbering(absId, relid);
1565  cellSM = relid[0];
1566 
1567  if (cellSM == sm) {
1568 
1569  if (fExcludeRejectedCells) {
1570  if (IsCellRejected(absId, kPHOS)) {
1571  continue;
1572  }
1573  }
1574 
1575  energy += phosCaloCells->GetCellAmplitude(absId);
1576  nCells += 1.;
1577 
1578  }
1579  }
1580  }
1581 
1582  if (returnNcells) {
1583  return nCells;
1584  }
1585 
1586  return energy;
1587 }
1588 
1593 {
1595  AliVCluster* cluster;
1596 
1597  UInt_t rejectionReason;
1598  Bool_t skipCell = kFALSE;
1599  Int_t clusType;
1600 
1602  for (AliClusterIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
1603 
1604  if (!clusCont->AcceptCluster(it.current_index(), rejectionReason)) {
1605 
1606  cluster = it->second;
1607 
1608  // check that the cell type and cluster type are the same
1609  if (cluster->IsEMCAL()) {
1610  clusType = kEMCal;
1611  }
1612  if (cluster->GetType() == AliVCluster::kPHOSNeutral) {
1613  clusType = kPHOS;
1614  }
1615  if (clusType != cellType) {
1616  continue;
1617  }
1618 
1619  // skip the cell if it belongs to a rejected cluster
1620  for (Int_t i = 0; i < cluster->GetNCells(); i++) {
1621 
1622  if (absId == cluster->GetCellAbsId(i)) {
1623  skipCell = kTRUE;
1624  }
1625 
1626  }
1627  }
1628  }
1629  return skipCell;
1630 }
1631 
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
const char * title
Definition: MakeQAPdf.C:27
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)
Int_t fNM02HistBins
! number of variable M02 bins
Double_t FindNearestNeighborDistance(AliTLorentzVector cluster)
THistManager fHistManager
Histogram manager.
AliPHOSGeometry * fPHOSGeo
! phos geometry
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.
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
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 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
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.
Bool_t fNeedEmcalGeom
whether or not the task needs the emcal geometry
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
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 coneType=0, Double_t R=0., Double_t Econe=0.)
static Double_t fgkEMCalDCalPhiDivide
phi value used to distinguish between DCal and EMCal