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