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