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