AliPhysics  80ccde44 (80ccde44)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskEmcalVsPhos.cxx
Go to the documentation of this file.
1 /**********************************************************************************
2 * Copyright (C) 2016, Copyright Holders of the ALICE Collaboration *
3 * All rights reserved. *
4 * *
5 * Redistribution and use in source and binary forms, with or without *
6 * modification, are permitted provided that the following conditions are met: *
7 * * Redistributions of source code must retain the above copyright *
8 * notice, this list of conditions and the following disclaimer. *
9 * * Redistributions in binary form must reproduce the above copyright *
10 * notice, this list of conditions and the following disclaimer in the *
11 * documentation and/or other materials provided with the distribution. *
12 * * Neither the name of the <organization> nor the *
13 * names of its contributors may be used to endorse or promote products *
14 * derived from this software without specific prior written permission. *
15 * *
16 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND *
17 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED *
18 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE *
19 * DISCLAIMED. IN NO EVENT SHALL ALICE COLLABORATION BE LIABLE FOR ANY *
20 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES *
21 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; *
22 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND *
23 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT *
24 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
25 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
26 * *********************************************************************************/
27 
28 #include <vector>
29 
30 #include <TClonesArray.h>
31 #include <TH1F.h>
32 #include <TH2F.h>
33 #include <TH3.h>
34 #include <TList.h>
35 #include <THnSparse.h>
36 
37 #include <AliVCluster.h>
38 #include <AliVParticle.h>
39 #include <AliLog.h>
40 
41 #include "AliTLorentzVector.h"
42 #include "AliEmcalJet.h"
43 #include "AliRhoParameter.h"
44 #include "AliJetContainer.h"
45 #include "AliParticleContainer.h"
46 #include "AliClusterContainer.h"
47 #include "AliEMCALGeometry.h"
48 #include "AliPHOSGeometry.h"
49 #include "AliOADBContainer.h"
50 
52 
56 
62  fHistManager(),
63  fEventCuts(0),
64  fEventCutList(0),
65  fUseManualEventCuts(kFALSE),
66  fUseAliEventCuts(kTRUE),
67  fMaxPt(200),
68  fNCentHistBins(0),
69  fCentHistBins(0),
70  fNPtHistBins(0),
71  fPtHistBins(0),
72  fNEtaBins(40),
73  fNPhiBins(200),
74  fPlotNeutralJets(kFALSE),
75  fPlotClustersInJets(kFALSE),
76  fPlotClusterHistograms(kTRUE),
77  fPlotCellHistograms(kTRUE),
78  fPlotClusWithoutNonLinCorr(kFALSE),
79  fPlotExotics(kFALSE),
80  fPHOSGeo(nullptr)
81 {
83 }
84 
91  AliAnalysisTaskEmcalJet(name, kTRUE),
92  fHistManager(name),
93  fEventCuts(0),
94  fEventCutList(0),
95  fUseManualEventCuts(kFALSE),
96  fUseAliEventCuts(kTRUE),
97  fMaxPt(200),
98  fNCentHistBins(0),
99  fCentHistBins(0),
100  fNPtHistBins(0),
101  fPtHistBins(0),
102  fNEtaBins(40),
103  fNPhiBins(200),
104  fPlotNeutralJets(kFALSE),
105  fPlotClustersInJets(kFALSE),
106  fPlotClusterHistograms(kTRUE),
107  fPlotCellHistograms(kTRUE),
108  fPlotClusWithoutNonLinCorr(kFALSE),
109  fPlotExotics(kFALSE),
110  fPHOSGeo(nullptr)
111 {
113 }
114 
119 {
120 }
121 
126 {
127  fNCentHistBins = 4;
129  fCentHistBins[0] = 0;
130  fCentHistBins[1] = 10;
131  fCentHistBins[2] = 30;
132  fCentHistBins[3] = 50;
133  fCentHistBins[4] = 90;
134 
135  fNPtHistBins = 82;
138  GenerateFixedBinArray(7, 0.3, 1, fPtHistBins+6);
139  GenerateFixedBinArray(10, 1, 3, fPtHistBins+13);
140  GenerateFixedBinArray(14, 3, 10, fPtHistBins+23);
141  GenerateFixedBinArray(10, 10, 20, fPtHistBins+37);
142  GenerateFixedBinArray(15, 20, 50, fPtHistBins+47);
143  GenerateFixedBinArray(20, 50, 150, fPtHistBins+62);
144 }
145 
151 {
153 
155 
156  TIter next(fHistManager.GetListOfHistograms());
157  TObject* obj = 0;
158  while ((obj = next())) {
159  fOutput->Add(obj);
160  }
161 
162  // Intialize AliEventCuts
163  if (fUseAliEventCuts) {
164  fEventCutList = new TList();
165  fEventCutList ->SetOwner();
166  fEventCutList ->SetName("EventCutOutput");
167 
168  fEventCuts.OverrideAutomaticTriggerSelection(fOffTrigger);
169  if(fUseManualEventCuts==1)
170  {
171  fEventCuts.SetManualMode();
172  // Configure manual settings here
173  // ...
174  }
175  fEventCuts.AddQAplotsToList(fEventCutList);
176  fOutput->Add(fEventCutList);
177  }
178 
179  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
180 }
181 
182 /*
183  * This function allocates the histograms for the calorimeter performance study.
184  */
186 {
187 
190  }
191 
192  if (fPlotCellHistograms) {
194  }
195 
196  if (fPlotNeutralJets) {
198  }
199 
200  if (fPlotClustersInJets) {
202  }
203 
204 }
205 
206 /*
207  * This function allocates the histograms for the calorimeter performance study.
208  */
210 {
211  TString histname;
212  TString htitle;
213  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
214 
215  Double_t* clusType = new Double_t[3+1];
216  GenerateFixedBinArray(3, -0.5, 2.5, clusType);
217  const Int_t nRejBins = 32;
218  Double_t* rejReasonBins = new Double_t[nRejBins+1];
219  GenerateFixedBinArray(nRejBins, 0, nRejBins, rejReasonBins);
220  const Int_t nExBins = 200;
221  Double_t* exBins = new Double_t[nExBins+1];
222  GenerateFixedBinArray(nExBins, 0, 1, exBins);
223 
224  AliEmcalContainer* cont = 0;
225  TIter nextClusColl(&fClusterCollArray);
226  while ((cont = static_cast<AliEmcalContainer*>(nextClusColl()))) {
227 
228  // rejection reason plot, to make efficiency correction
229  histname = TString::Format("%s/hClusterRejectionReasonEMCal", cont->GetArrayName().Data());
230  htitle = histname + ";Rejection reason;#it{E}_{clus} (GeV/)";
231  TH2* hist = fHistManager.CreateTH2(histname.Data(), htitle.Data(), nRejBins, rejReasonBins, fNPtHistBins, fPtHistBins);
232  SetRejectionReasonLabels(hist->GetXaxis());
233 
234  histname = TString::Format("%s/hClusterRejectionReasonPHOS", cont->GetArrayName().Data());
235  htitle = histname + ";Rejection reason;#it{E}_{clus} (GeV/)";
236  TH2* histPhos = fHistManager.CreateTH2(histname.Data(), htitle.Data(), nRejBins, rejReasonBins, fNPtHistBins, fPtHistBins);
237  SetRejectionReasonLabels(histPhos->GetXaxis());
238 
239  // plot by SM
240  const Int_t nEmcalSM = 20;
241  for (Int_t sm = 0; sm < nEmcalSM; sm++) {
242  histname = TString::Format("%s/BySM/hEmcalClusEnergy_SM%d", cont->GetArrayName().Data(), sm);
243  htitle = histname + ";#it{E}_{cluster} (GeV);counts";
244  fHistManager.CreateTH1(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins);
245  }
246 
247  for (Int_t sm = 1; sm < 5; sm++) {
248  histname = TString::Format("%s/BySM/hPhosClusEnergy_SM%d", cont->GetArrayName().Data(), sm);
249  htitle = histname + ";#it{E}_{cluster} (GeV);counts";
250  fHistManager.CreateTH1(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins);
251  }
252 
253  // Plot cluster THnSparse (centrality, cluster type, E, E-hadcorr, has matched track, M02, Ncells)
254  Int_t dim = 0;
255  TString title[20];
256  Int_t nbins[20] = {0};
257  Double_t min[30] = {0.};
258  Double_t max[30] = {0.};
259  Double_t *binEdges[20] = {0};
260 
262  title[dim] = "Centrality %";
263  nbins[dim] = fNCentHistBins;
264  binEdges[dim] = fCentHistBins;
265  min[dim] = fCentHistBins[0];
266  max[dim] = fCentHistBins[fNCentHistBins];
267  dim++;
268  }
269 
270  title[dim] = "#eta";
271  nbins[dim] = 28;
272  min[dim] = -0.7;
273  max[dim] = 0.7;
274  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
275  dim++;
276 
277  title[dim] = "#phi";
278  nbins[dim] = 100;
279  min[dim] = 1.;
280  max[dim] = 6.;
281  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
282  dim++;
283 
284  title[dim] = "#it{E}_{clus} (GeV)";
285  nbins[dim] = fNPtHistBins;
286  binEdges[dim] = fPtHistBins;
287  min[dim] = fPtHistBins[0];
288  max[dim] = fPtHistBins[fNPtHistBins];
289  dim++;
290 
291  title[dim] = "#it{E}_{clus, hadcorr} or #it{E}_{core} (GeV)";
292  nbins[dim] = fNPtHistBins;
293  binEdges[dim] = fPtHistBins;
294  min[dim] = fPtHistBins[0];
295  max[dim] = fPtHistBins[fNPtHistBins];
296  dim++;
297 
298  title[dim] = "Matched track";
299  nbins[dim] = 2;
300  min[dim] = -0.5;
301  max[dim] = 1.5;
302  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
303  dim++;
304 
305  title[dim] = "M02";
306  nbins[dim] = 50;
307  min[dim] = 0;
308  max[dim] = 5;
309  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
310  dim++;
311 
312  title[dim] = "Ncells";
313  nbins[dim] = 30;
314  min[dim] = -0.5;
315  max[dim] = 29.5;
316  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
317  dim++;
318 
319  title[dim] = "Dispersion cut";
320  nbins[dim] = 2;
321  min[dim] = -0.5;
322  max[dim] = 1.5;
323  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
324  dim++;
325 
326  TString thnname = TString::Format("%s/clusterObservables", cont->GetArrayName().Data());
327  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
328  for (Int_t i = 0; i < dim; i++) {
329  hn->GetAxis(i)->SetTitle(title[i]);
330  hn->SetBinEdges(i, binEdges[i]);
331  }
332 
333  if (fPlotExotics) {
334  histname = TString::Format("%s/hFcrossEMCal", cont->GetArrayName().Data());
335  htitle = histname + ";Centrality (%);Fcross;#it{E}_{clus} (GeV/)";
336  TH3* hist = fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, nExBins, exBins, fNPtHistBins, fPtHistBins);
337  }
338  }
339 
340 }
341 
342 /*
343  * This function allocates the histograms for calo cells.
344  */
346 {
347  TString histname;
348  TString htitle;
349 
350  Double_t* clusType = new Double_t[3+1];
351  GenerateFixedBinArray(3, -0.5, 2.5, clusType);
352 
353  // centrality vs. cell energy vs. cell type (for all cells)
354  histname = TString::Format("Cells/hCellEnergyAll");
355  htitle = histname + ";#it{E}_{cell} (GeV);Centrality (%); Cluster type";
356  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins, 3, clusType);
357 
358  // centrality vs. cell energy vs. cell type (for cells in accepted clusters)
359  histname = TString::Format("Cells/hCellEnergyAccepted");
360  htitle = histname + ";#it{E}_{cell} (GeV);Centrality (%); Cluster type";
361  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins, 3, clusType);
362 
363  // centrality vs. cell energy vs. cell type (for leading cells in accepted clusters)
364  histname = TString::Format("Cells/hCellEnergyLeading");
365  htitle = histname + ";#it{E}_{cell} (GeV);Centrality (%); Cluster type";
366  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins, 3, clusType);
367 
368  // plot cell patches by SM
369  const Int_t nEmcalSM = 20;
370  for (Int_t sm = 0; sm < nEmcalSM; sm++) {
371  histname = TString::Format("Cells/BySM/hEmcalPatchEnergy_SM%d", sm);
372  htitle = histname + ";#it{E}_{cell patch} (GeV);Centrality (%)";
373  fHistManager.CreateTH2(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins);
374  }
375 
376  for (Int_t sm = 1; sm < 5; sm++) {
377  histname = TString::Format("Cells/BySM/hPhosPatchEnergy_SM%d", sm);
378  htitle = histname + ";#it{E}_{cell patch} (GeV);Centrality (%)";
379  fHistManager.CreateTH2(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins);
380  }
381 
382 }
383 
384 /*
385  * This function allocates the histograms for neutral jets.
386  */
388 {
389  AliJetContainer* jets = 0;
390  TIter nextJetColl(&fJetCollArray);
391  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
392 
393  TString axisTitle[30]= {""};
394  Int_t nbinsJet[30] = {0};
395  Double_t minJet[30] = {0.};
396  Double_t maxJet[30] = {0.};
397  Double_t *binEdgesJet[20] = {0};
398  Int_t dimJet = 0;
399 
400  if (fForceBeamType != kpp) {
401  axisTitle[dimJet] = "Centrality (%)";
402  nbinsJet[dimJet] = fNCentHistBins;
403  binEdgesJet[dimJet] = fCentHistBins;
404  minJet[dimJet] = fCentHistBins[0];
405  maxJet[dimJet] = fCentHistBins[fNCentHistBins];
406  dimJet++;
407  }
408 
409  axisTitle[dimJet] = "#eta_{jet}";
410  nbinsJet[dimJet] = 28;
411  minJet[dimJet] = -0.7;
412  maxJet[dimJet] = 0.7;
413  binEdgesJet[dimJet] = GenerateFixedBinArray(nbinsJet[dimJet], minJet[dimJet], maxJet[dimJet]);
414  dimJet++;
415 
416  axisTitle[dimJet] = "#phi_{jet} (rad)";
417  nbinsJet[dimJet] = 100;
418  minJet[dimJet] = 1.;
419  maxJet[dimJet] = 6.;
420  binEdgesJet[dimJet] = GenerateFixedBinArray(nbinsJet[dimJet], minJet[dimJet], maxJet[dimJet]);
421  dimJet++;
422 
423  axisTitle[dimJet] = "#it{E}_{T} (GeV)";
424  nbinsJet[dimJet] = fNPtHistBins;
425  binEdgesJet[dimJet] = fPtHistBins;
426  minJet[dimJet] = fPtHistBins[0];
427  maxJet[dimJet] = fPtHistBins[fNPtHistBins];
428  dimJet++;
429 
430  axisTitle[dimJet] = "#rho (GeV/#it{c})";
431  nbinsJet[dimJet] = 100;
432  minJet[dimJet] = 0.;
433  maxJet[dimJet] = 1000.;
434  binEdgesJet[dimJet] = GenerateFixedBinArray(nbinsJet[dimJet], minJet[dimJet], maxJet[dimJet]);
435  dimJet++;
436 
437  axisTitle[dimJet] = "N_{clusters}";
438  nbinsJet[dimJet] = 20;
439  minJet[dimJet] = -0.5;
440  maxJet[dimJet] = 19.5;
441  binEdgesJet[dimJet] = GenerateFixedBinArray(nbinsJet[dimJet], minJet[dimJet], maxJet[dimJet]);
442  dimJet++;
443 
444  TString thnname = TString::Format("%s/hClusterJetObservables", jets->GetArrayName().Data());
445  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dimJet, nbinsJet, minJet, maxJet);
446  for (Int_t i = 0; i < dimJet; i++) {
447  hn->GetAxis(i)->SetTitle(axisTitle[i]);
448  hn->SetBinEdges(i, binEdgesJet[i]);
449  }
450 
451  }
452 }
453 
454 /*
455  * This function allocates the histograms for clusters within jets.
456  */
458 {
459  TString histname;
460  TString htitle;
461  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
462 
463  AliJetContainer* jets = 0;
464  TIter nextJetColl(&fJetCollArray);
465  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
466 
467  // Plot cluster spectra within jets
468  // (centrality, cluster energy, jet pT, jet eta, jet phi)
469  Int_t dim = 0;
470  TString title[20];
471  Int_t nbins[20] = {0};
472  Double_t min[30] = {0.};
473  Double_t max[30] = {0.};
474  Double_t *binEdges[20] = {0};
475 
477  title[dim] = "Centrality %";
478  nbins[dim] = fNCentHistBins;
479  binEdges[dim] = fCentHistBins;
480  min[dim] = fCentHistBins[0];
481  max[dim] = fCentHistBins[fNCentHistBins];
482  dim++;
483  }
484 
485  title[dim] = "#it{E}_{clus} (GeV)";
486  nbins[dim] = fNPtHistBins;
487  binEdges[dim] = fPtHistBins;
488  min[dim] = fPtHistBins[0];
489  max[dim] = fPtHistBins[fNPtHistBins];
490  dim++;
491 
492  title[dim] = "#it{p}_{T,jet}^{corr}";
493  nbins[dim] = nPtBins;
494  min[dim] = 0;
495  max[dim] = fMaxPt;
496  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
497  dim++;
498 
499  title[dim] = "#eta_{jet}";
500  nbins[dim] = fNEtaBins;
501  min[dim] = -0.5;
502  max[dim] = 0.5;
503  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
504  dim++;
505 
506  title[dim] = "#phi_{jet}";
507  nbins[dim] = fNPhiBins;
508  min[dim] = 1.;
509  max[dim] = 6.;
510  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
511  dim++;
512 
513  TString thnname = TString::Format("%s/hClustersInJets", jets->GetArrayName().Data());
514  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
515  for (Int_t i = 0; i < dim; i++) {
516  hn->GetAxis(i)->SetTitle(title[i]);
517  hn->SetBinEdges(i, binEdges[i]);
518  }
519 
520  // (jet type, jet pT, cluster shift)
521  histname = TString::Format("%s/hCaloJESshift", jets->GetArrayName().Data());
522  htitle = histname + ";type;#it{p}_{T}^{corr} (GeV/#it{c});#Delta_{JES}";
523  fHistManager.CreateTH3(histname.Data(), htitle.Data(), 3, -0.5, 2.5, nPtBins, 0, fMaxPt, 100, 0, 20);
524 
525  }
526 }
527 
533 {
535 
536  fNeedEmcalGeom = kTRUE;
537 
538  // Load the PHOS geometry
539  fPHOSGeo = AliPHOSGeometry::GetInstance();
540  if (fPHOSGeo) {
541  AliInfo("Found instance of PHOS geometry!");
542  }
543  else {
544  AliInfo("Creating PHOS geometry!");
545  Int_t runNum = InputEvent()->GetRunNumber();
546  if(runNum<209122) //Run1
547  fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP");
548  else
549  fPHOSGeo = AliPHOSGeometry::GetInstance("Run2");
550 
551  if (fPHOSGeo) {
552  AliOADBContainer geomContainer("phosGeo");
553  geomContainer.InitFromFile("$ALICE_PHYSICS/OADB/PHOS/PHOSMCGeometry.root","PHOSMCRotationMatrixes");
554  TObjArray* matrixes = (TObjArray*)geomContainer.GetObject(runNum,"PHOSRotationMatrixes");
555  for(Int_t mod=0; mod<6; mod++) {
556  if(!matrixes->At(mod)) continue;
557  fPHOSGeo->SetMisalMatrix(((TGeoHMatrix*)matrixes->At(mod)),mod);
558  printf(".........Adding Matrix(%d), geo=%p\n",mod,fPHOSGeo);
559  ((TGeoHMatrix*)matrixes->At(mod))->Print();
560  }
561  }
562  }
563 }
564 
569 {
570  if (fUseAliEventCuts) {
571  if (!fEventCuts.AcceptEvent(InputEvent()))
572  {
573  PostData(1, fOutput);
574  return kFALSE;
575  }
576  }
577  else {
579  }
580  return kTRUE;
581 }
582 
591 {
592  return kTRUE;
593 }
594 
602 {
603 
606  }
607 
608  if (fPlotNeutralJets) {
610  }
611 
612  if (fPlotClustersInJets) {
614  }
615 
616  if (fPlotCellHistograms) {
618  }
619 
620  return kTRUE;
621 }
622 
623 /*
624  * This function fills the histograms for the calorimeter performance study.
625  */
627 {
628  // Define some vars
629  TString histname;
630  Double_t Enonlin;
631  Double_t Ehadcorr;
632  Int_t absId;
633  Double_t ecell;
634  Double_t leadEcell;
635 
636  // Get cells from event
637  fCaloCells = InputEvent()->GetEMCALCells();
638  AliVCaloCells* phosCaloCells = InputEvent()->GetPHOSCells();
639 
640  // Loop through clusters and plot cluster THnSparse (centrality, cluster type, E, E-hadcorr, has matched track, M02, Ncells)
641  AliClusterContainer* clusters = 0;
642  TIter nextClusColl(&fClusterCollArray);
643  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
645  for (AliClusterIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
646 
647  // Determine cluster type (EMCal/DCal/Phos)
648  ClusterType clusType = kNA;
649  if (it->second->IsEMCAL()) {
650  Double_t phi = it->first.Phi_0_2pi();
651  Int_t isDcal = Int_t(phi > fgkEMCalDCalPhiDivide);
652  if (isDcal == 0) {
653  clusType = kEMCal;
654  } else if (isDcal == 1) {
655  clusType = kDCal;
656  }
657  } else if (it->second->GetType() == AliVCluster::kPHOSNeutral){
658  clusType = kPHOS;
659  }
660 
661  // rejection reason plots, to make efficiency correction
662  if (it->second->IsEMCAL()) {
663  histname = TString::Format("%s/hClusterRejectionReasonEMCal", clusters->GetArrayName().Data());
664  UInt_t rejectionReason = 0;
665  if (!clusters->AcceptCluster(it.current_index(), rejectionReason)) {
666  fHistManager.FillTH2(histname, clusters->GetRejectionReasonBitPosition(rejectionReason), it->first.E());
667  continue;
668  }
669  } else if (it->second->GetType() == AliVCluster::kPHOSNeutral){
670  histname = TString::Format("%s/hClusterRejectionReasonPHOS", clusters->GetArrayName().Data());
671  UInt_t rejectionReason = 0;
672  if (!clusters->AcceptCluster(it.current_index(), rejectionReason)) {
673  fHistManager.FillTH2(histname, clusters->GetRejectionReasonBitPosition(rejectionReason), it->first.E());
674  continue;
675  }
676  } else {
677  continue;
678  }
679 
680  // Fill cluster spectra by SM, and fill cell histograms
681  Enonlin = 0;
682  Ehadcorr = 0;
683  if (it->second->IsEMCAL()) {
684 
685  Ehadcorr = it->second->GetHadCorrEnergy();
686  Enonlin = it->second->GetNonLinCorrEnergy();
688  Enonlin = it->second->E();
689  }
690 
691  if (fPlotExotics) {
692  histname = TString::Format("%s/hFcrossEMCal", clusters->GetArrayName().Data());
693  Double_t Fcross = GetFcross(it->second, fCaloCells);
694  fHistManager.FillTH3(histname, fCent, Fcross, it->second->E());
695  }
696 
697  Int_t sm = fGeom->GetSuperModuleNumber(it->second->GetCellAbsId(0));
698  if (sm >=0 && sm < 20) {
699  histname = TString::Format("%s/BySM/hEmcalClusEnergy_SM%d", clusters->GetArrayName().Data(), sm);
700  fHistManager.FillTH1(histname, it->second->E());
701  }
702  else {
703  AliError(Form("Supermodule %d does not exist!", sm));
704  }
705 
706  // Get cells from each accepted cluster, and plot centrality vs. cell energy vs. cell type
707  histname = TString::Format("Cells/hCellEnergyAccepted");
708  leadEcell = 0;
709  for (Int_t iCell = 0; iCell < it->second->GetNCells(); iCell++){
710  absId = it->second->GetCellAbsId(iCell);
711  ecell = fCaloCells->GetCellAmplitude(absId);
712  fHistManager.FillTH3(histname, ecell, fCent, kEMCal); // Note: I don't distinguish EMCal from DCal cells
713  if (ecell > leadEcell) {
714  leadEcell = ecell;
715  }
716  }
717  // Plot also the leading cell
718  histname = TString::Format("Cells/hCellEnergyLeading");
719  fHistManager.FillTH3(histname, leadEcell, fCent, kEMCal);
720 
721  } else if (it->second->GetType() == AliVCluster::kPHOSNeutral){
722 
723  Ehadcorr = it->second->GetCoreEnergy();
724  Enonlin = it->second->E();
725 
726  Int_t relid[4];
727  if (fPHOSGeo) {
728  fPHOSGeo->AbsToRelNumbering(it->second->GetCellAbsId(0), relid);
729  Int_t sm = relid[0];
730  if (sm >=1 && sm < 5) {
731  histname = TString::Format("%s/BySM/hPhosClusEnergy_SM%d", clusters->GetArrayName().Data(), sm);
732  fHistManager.FillTH1(histname, it->second->E());
733  }
734  else {
735  AliError(Form("Supermodule %d does not exist!", sm));
736  }
737  }
738 
739  // Get cells from each accepted cluster, and plot centrality vs. cell energy vs. cell type
740  histname = TString::Format("Cells/hCellEnergyAccepted");
741  leadEcell = 0;
742  for (Int_t iCell = 0; iCell < it->second->GetNCells(); iCell++){
743  absId = it->second->GetCellAbsId(iCell);
744  ecell = phosCaloCells->GetCellAmplitude(absId);
745  fHistManager.FillTH3(histname, ecell, fCent, kPHOS);
746  if (ecell > leadEcell) {
747  leadEcell = ecell;
748  }
749  }
750  // Plot also the leading cell
751  histname = TString::Format("Cells/hCellEnergyLeading");
752  fHistManager.FillTH3(histname, leadEcell, fCent, kPHOS);
753  }
754 
755  // Check if the cluster has a matched track
756  Int_t hasMatchedTrack = -1;
757  Int_t nMatchedTracks = it->second->GetNTracksMatched();
758  if (nMatchedTracks == 0) {
759  hasMatchedTrack = 0;
760  } else if (nMatchedTracks > 0) {
761  hasMatchedTrack = 1;
762  }
763 
764  // Check if the cluster passes the dispersion cut for photon-like cluster (meaningful only for PHOS)
765  Int_t passedDispersionCut = 0;
766  if (it->second->Chi2() < 2.5*2.5) {
767  passedDispersionCut = 1;
768  }
769 
770  Double_t contents[30]={0};
771  histname = TString::Format("%s/clusterObservables", clusters->GetArrayName().Data());
772  THnSparse* histClusterObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
773  if (!histClusterObservables) return;
774  for (Int_t i = 0; i < histClusterObservables->GetNdimensions(); i++) {
775  TString title(histClusterObservables->GetAxis(i)->GetTitle());
776  if (title=="Centrality %")
777  contents[i] = fCent;
778  else if (title=="#eta")
779  contents[i] = it->first.Eta();
780  else if (title=="#phi")
781  contents[i] = it->first.Phi_0_2pi();
782  else if (title=="#it{E}_{clus} (GeV)")
783  contents[i] = Enonlin;
784  else if (title=="#it{E}_{clus, hadcorr} or #it{E}_{core} (GeV)")
785  contents[i] = Ehadcorr;
786  else if (title=="Matched track")
787  contents[i] = hasMatchedTrack;
788  else if (title=="M02")
789  contents[i] = it->second->GetM02();
790  else if (title=="Ncells")
791  contents[i] = it->second->GetNCells();
792  else if (title=="Dispersion cut")
793  contents[i] = passedDispersionCut;
794  else
795  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
796  }
797  histClusterObservables->Fill(contents);
798 
799  }
800 
801  }
802 }
803 
804 /*
805  * This function fills cell histograms.
806  */
808 {
809  // Define some vars
810  TString histname;
811  Int_t absId;
812  Double_t ecell;
813 
814  // Get cells from event
815  fCaloCells = InputEvent()->GetEMCALCells();
816  AliVCaloCells* phosCaloCells = InputEvent()->GetPHOSCells();
817 
818  // Loop through all cells and fill histos
819  Int_t sm;
820  Int_t relid[4];
821  Double_t patchSumEMCal[20] = {0.};
822  Double_t patchSumPHOS[4] = {0.};
823  for (Int_t i=0; i<fCaloCells->GetNumberOfCells(); i++) {
824 
825  absId = fCaloCells->GetCellNumber(i);
826  ecell = fCaloCells->GetCellAmplitude(absId);
827 
828  // Fill cell histo
829  histname = TString::Format("Cells/hCellEnergyAll");
830  fHistManager.FillTH3(histname, ecell, fCent, kEMCal); // Note: I don't distinguish EMCal from DCal cells
831 
832  // Fill cell patch histo, per SM
833  sm = fGeom->GetSuperModuleNumber(absId);
834  if (sm >=0 && sm < 20) {
835  patchSumEMCal[sm] += ecell;
836  }
837 
838  }
839 
840  for (Int_t i=0; i<phosCaloCells->GetNumberOfCells(); i++) {
841 
842  absId = phosCaloCells->GetCellNumber(i);
843  ecell = phosCaloCells->GetCellAmplitude(absId);
844 
845  // Fill cell histo
846  histname = TString::Format("Cells/hCellEnergyAll");
847  fHistManager.FillTH3(histname, ecell, fCent, kPHOS);
848 
849  // Fill cell patch histo, per SM
850  fPHOSGeo->AbsToRelNumbering(absId, relid);
851  sm = relid[0];
852  if (sm >=1 && sm < 5) {
853  patchSumPHOS[sm-1] += ecell;
854  }
855 
856  }
857 
858  for (Int_t sm = 0; sm < 20; sm++) {
859  histname = TString::Format("Cells/BySM/hEmcalPatchEnergy_SM%d", sm);
860  fHistManager.FillTH2(histname, patchSumEMCal[sm], fCent);
861  }
862 
863  for (Int_t sm = 1; sm < 5; sm++) {
864  histname = TString::Format("Cells/BySM/hPhosPatchEnergy_SM%d", sm);
865  fHistManager.FillTH2(histname, patchSumPHOS[sm-1], fCent);
866  }
867 
868 }
869 
870 /*
871  * This function fills neutral jet histograms.
872  */
874 {
875  AliJetContainer* jets = 0;
876  TIter nextJetColl(&fJetCollArray);
877  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
878 
879  // plot neutral jets THnSparse (centrality, eta, phi, E, Nclusters)
880  Double_t contents[30]={0};
881  TString histname = TString::Format("%s/hClusterJetObservables", jets->GetArrayName().Data());
882  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
883  if (!histJetObservables) return;
884 
885  for (auto jet : jets->accepted()) {
886 
887  for (Int_t i = 0; i < histJetObservables->GetNdimensions(); i++) {
888  TString title(histJetObservables->GetAxis(i)->GetTitle());
889  if (title=="Centrality (%)")
890  contents[i] = fCent;
891  else if (title=="#eta_{jet}")
892  contents[i] = jet->Eta();
893  else if (title=="#phi_{jet} (rad)")
894  contents[i] = jet->Phi_0_2pi();
895  else if (title=="#it{E}_{T} (GeV)")
896  contents[i] = jet->Pt();
897  else if (title=="#rho (GeV/#it{c})")
898  contents[i] = jet->Pt() / jet->Area();
899  else if (title=="N_{clusters}")
900  contents[i] = jet->GetNumberOfClusters();
901  else
902  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
903  }
904  histJetObservables->Fill(contents);
905 
906  }
907  }
908 }
909 
910 /*
911  * This function fills clusters within jets and estimates the JES shift due to the bump.
912  */
914 {
915  TString histname;
916  AliJetContainer* jets = 0;
917  TIter nextJetColl(&fJetCollArray);
918  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
919 
920  for (auto jet : jets->accepted()) {
921 
922  // Fill cluster spectra of clusters within jets
923  //(centrality, cluster energy, jet pT, jet eta, jet phi)
924  histname = TString::Format("%s/hClustersInJets", jets->GetArrayName().Data());
925  Int_t nClusters = jet->GetNumberOfClusters();
926  AliVCluster* clus;
927  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
928  clus = jet->Cluster(iClus);
929  Double_t x[5] = {fCent, clus->E(), GetJetPt(jets, jet), jet->Eta(), jet->Phi_0_2pi()};
930  fHistManager.FillTHnSparse(histname, x);
931  }
932 
933  // Loop through clusters, and plot estimated shift in JES due to cluster bump
934  // Only do for 0-10% centrality, and for EMCal/DCal
935  Double_t eclus;
936  Double_t shift;
937  Double_t shiftSum = 0;
938  if (fCent < 10.) {
939  if (GetJetType(jet) > -0.5 && GetJetType(jet) < 1.5) {
940  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
941  clus = jet->Cluster(iClus);
942  eclus = clus->E();
943  if (eclus > 0.5) {
944  shift = 0.79 * TMath::Exp(-0.5 * ((eclus - 3.81) / 1.50)*((eclus - 3.81) / 1.50) );
945  shiftSum += shift;
946  }
947  }
948  histname = TString::Format("%s/hCaloJESshift", jets->GetArrayName().Data());
949  fHistManager.FillTH3(histname, GetJetType(jet), GetJetPt(jets, jet), shiftSum);
950  }
951  }
952 
953  }
954  }
955 }
956 
961 {
962  // Get eta-phi dependent jet pT scale factor
963  Double_t jetPt = jet->Pt();
964 
965  // Compute pTcorr
966  Double_t rho = jetCont->GetRhoVal();
967  Double_t pT = jetPt - rho * jet->Area();
968 
969  return pT;
970 }
971 
976 {
977  Double_t deltaPhi = TMath::Abs(part->Phi_0_2pi() - phiRef);
978  Double_t deltaEta = TMath::Abs(part->Eta() - etaRef);
979  Double_t deltaR = TMath::Sqrt( deltaPhi*deltaPhi + deltaEta*deltaEta );
980  return deltaR;
981 }
982 
987 {
988  UInt_t jetType = jet->GetJetAcceptanceType();
989  Double_t type = -1;
990  if (jetType & AliEmcalJet::kEMCAL) {
991  type = kEMCal;
992  }
993  else if (jetType & AliEmcalJet::kDCALonly) {
994  type = kDCal;
995  }
996  else if (jetType & AliEmcalJet::kPHOS) {
997  type = kPHOS;
998  }
999 
1000  return type;
1001 }
1002 
1006 Double_t AliAnalysisTaskEmcalVsPhos::GetFcross(AliVCluster *cluster, AliVCaloCells *cells)
1007 {
1008  Int_t AbsIdseed = -1;
1009  Double_t Eseed = 0;
1010  for (Int_t i = 0; i < cluster->GetNCells(); i++) {
1011  if (cells->GetCellAmplitude(cluster->GetCellAbsId(i)) > Eseed) {
1012  Eseed = cells->GetCellAmplitude(cluster->GetCellAbsId(i));
1013  AbsIdseed = cluster->GetCellAbsId(i);
1014  }
1015  }
1016 
1017  if (Eseed < 1e-9) {
1018  return 100;
1019  }
1020 
1021  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
1022  fGeom->GetCellIndex(AbsIdseed,imod,iTower,iIphi,iIeta);
1023  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi,iIeta,iphi,ieta);
1024 
1025  //Get close cells index and energy, not in corners
1026 
1027  Int_t absID1 = -1;
1028  Int_t absID2 = -1;
1029 
1030  if (iphi < AliEMCALGeoParams::fgkEMCALRows-1) {
1031  absID1 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
1032  }
1033  if (iphi > 0) {
1034  absID2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
1035  }
1036 
1037  // In case of cell in eta = 0 border, depending on SM shift the cross cell index
1038 
1039  Int_t absID3 = -1;
1040  Int_t absID4 = -1;
1041 
1042  if (ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2)) {
1043  absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
1044  absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
1045  }
1046  else if (ieta == 0 && imod%2) {
1047  absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
1048  absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
1049  }
1050  else {
1051  if (ieta < AliEMCALGeoParams::fgkEMCALCols-1) {
1052  absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
1053  }
1054  if (ieta > 0) {
1055  absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
1056  }
1057  }
1058 
1059  Double_t ecell1 = cells->GetCellAmplitude(absID1);
1060  Double_t ecell2 = cells->GetCellAmplitude(absID2);
1061  Double_t ecell3 = cells->GetCellAmplitude(absID3);
1062  Double_t ecell4 = cells->GetCellAmplitude(absID4);
1063 
1064  Double_t Ecross = ecell1 + ecell2 + ecell3 + ecell4;
1065 
1066  Double_t Fcross = 1 - Ecross/Eseed;
1067 
1068  return Fcross;
1069 }
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
Double_t GetJetPt(AliJetContainer *jetCont, AliEmcalJet *jet)
const char * title
Definition: MakeQAPdf.C:27
Definition: External.C:244
PHOS acceptance.
Definition: AliEmcalJet.h:68
Bool_t fPlotClustersInJets
Set whether to plot histogram of clusters within jets.
UInt_t fOffTrigger
offline trigger for event selection
Int_t fNPhiBins
Number of phi bins in DCal region (for background/correction)
bidirectional stl iterator over the EMCAL iterable container
Int_t fNCentHistBins
! number of cent bins
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)
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="")
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
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
Double_t Phi_0_2pi() const
EMCal acceptance.
Definition: AliEmcalJet.h:62
Double_t GetJetType(AliEmcalJet *jet)
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
BeamType fForceBeamType
forced beam type
virtual Bool_t AcceptCluster(Int_t i, UInt_t &rejectionReason) const
Double_t fCent
!event centrality
Double_t GetFcross(AliVCluster *cluster, AliVCaloCells *cells)
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
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
Double_t GetDeltaR(AliTLorentzVector *part, Double_t etaRef, Double_t phiRef)
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)
AliEmcalList * fOutput
!output list
Definition: External.C:220
Int_t fNEtaBins
Number of eta bins in DCal region (for background/correction)
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
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
static Double_t fgkEMCalDCalPhiDivide
phi value used to distinguish between DCal and EMCal