AliPhysics  93c564e (93c564e)
AliAnalysisTaskEmcalRecalcPatchesRef.cxx
Go to the documentation of this file.
1 /************************************************************************************
2  * Copyright (C) 2017, 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 #include <algorithm>
28 #include <array>
29 #include <iostream>
30 #include <unordered_map>
31 #include <sstream>
32 #include <string>
33 #include <TClonesArray.h>
34 #include <THistManager.h>
35 #include <TLinearBinning.h>
36 #include <TList.h>
37 #include <TString.h>
38 #include "AliAnalysisManager.h"
39 #include "AliEMCALTriggerPatchInfo.h"
42 
44 
45 using namespace EMCalTriggerPtAnalysis;
46 
49  fEnableSumw2(false),
50  fOnlineThresholds()
51 {
52 
53 }
54 
57  fEnableSumw2(false),
59 {
60  fOnlineThresholds.Set(8);
61 }
62 
64  const std::array<std::string, 9> kNamesTriggerClasses = {{"MB", "EG1", "EG2", "DG1", "DG2", "EJ1", "EJ2", "DJ1", "DJ2"}};
65  const std::array<std::string, 6> kNamesTriggerClusters = {{"ANY", "CENT", "CENTNOTRD", "BOTH", "ONLYCENT", "ONLYCENTNOTRD"}};
66  const std::array<std::string, 4> kNamesPatchTypes = {{"EGA", "DGA", "EJE", "DJE"}};
67 
68  TLinearBinning ADCBinning(2000, 0., 2000.), colbinning(48, -0.5, 47.5), rowbinning(104, -0.5, 103.5), npatchbinning(51, -0.5, 50.5), noverlapbinning(21, -0.5, 20.5);
69  const TBinning *firedpatchbinning[5] = {&ADCBinning, &colbinning, &rowbinning, &npatchbinning, &noverlapbinning},
70  *allpatchbinning[3] = {&ADCBinning, &npatchbinning, &noverlapbinning};
71  // Create event counters
72  for(const auto &kt : kNamesTriggerClasses){
73  // Min. Bias: Only CENT cluster - no distinction between trigger clusters
74  // EMCAL triggers: Distinction between CENT and CENTNOTRD necessary
75  if(kt.find("MB") != std::string::npos)
76  fHistos->CreateTH1(Form("hEventCounter%s", kt.data()), Form("Event counter for %s", kt.data()), 1, 0.5, 1.5);
77  else {
78  for(const auto &kc : kNamesTriggerClusters) {
79  fHistos->CreateTH1(Form("hEventCounter%s%s", kt.data(), kc.data()), Form("Event counter for %s in cluster %s", kt.data(), kc.data()), 1, 0.5, 1.5);
80  }
81  }
82  }
83 
84  // Min. Bias: Create patch energy spectra (all patches) for EMCAL and DCAL
85  // Min. Bias trigger only in CENT cluster
86  for(const auto &kp : kNamesPatchTypes) fHistos->CreateTH1(Form("hPatchADC%sMB", kp.data()), Form("Patch ADC spectra for %s patches in MB events", kp.data()), 2000, 0., 2000., fEnableSumw2 ? "s" : "");
87 
88  // Triggers: Create trigger spectra and THnSparse of firing patches
89  for(const auto &kt : kNamesTriggerClasses) {
90  if(kt == "MB") continue;
91  const char detector = kt[0];
92  const char *patchtype = ((kt[1] == 'G') ? "GA" : "JE");
93  // distinction between trigger clusters
94  for(const auto &kc : kNamesTriggerClusters){
95  fHistos->CreateTH1(Form("hPatchADC%c%s%s%s", detector, patchtype, kt.data(), kc.data()), Form("Patch ADC spectra for %c%s patches in %s events (cluster %s)", detector, patchtype, kt.data(), kc.data()), 2000, 0., 2000., fEnableSumw2 ? "s" : "");
96  fHistos->CreateTHnSparse(Form("hFiredPatches%c%s%s%s", detector, patchtype, kt.data(), kc.data()), Form("Fired %c%s patches for trigger %s (cluster %s)", detector, patchtype, kt.data(), kc.data()), 5, firedpatchbinning, fEnableSumw2 ? "s" : "");
97  fHistos->CreateTHnSparse(Form("hAllPatches%c%s%s%s", detector, patchtype, kt.data(), kc.data()), Form("Fired %c%s patches for trigger %s (cluster %s)", detector, patchtype, kt.data(), kc.data()), 3, allpatchbinning, fEnableSumw2 ? "s" : "");
98  }
99  }
100 }
101 
103  const std::array<std::string, 9> kNamesTriggerClasses = {{"MB", "EG1", "EG2", "DG1", "DG2", "EJ1", "EJ2", "DJ1", "DJ2"}};
104  const auto selclusters = GetAcceptedTriggerClusters(fInputEvent->GetFiredTriggerClasses().Data());
105  for(const auto &kt : kNamesTriggerClasses) {
106  if(std::find(fSelectedTriggers.begin(), fSelectedTriggers.end(), kt) != fSelectedTriggers.end()) {
107  if(kt == "MB") fHistos->FillTH1(Form("hEventCounter%s", kt.data()), 1.);
108  else {
109  for(const auto &kc : selclusters) fHistos->FillTH1(Form("hEventCounter%s%s", kt.data(), kc.data()), 1.);
110  }
111  }
112  }
113 }
114 
116  const std::array<std::string, 9> kNamesTriggerClasses = {{"MB", "EG1", "EG2", "DG1", "DG2", "EJ1", "EJ2", "DJ1", "DJ2"}};
117  const std::unordered_map<std::string, ETriggerThreshold_t> kPatchIndex = {{"EG1", kThresholdEG1},{"EG2", kThresholdEG2},{"DG1", kThresholdDG1},{"DG2", kThresholdDG2},
118  {"EJ1", kThresholdEJ1},{"EJ2", kThresholdEJ2},{"DJ1", kThresholdDJ1},{"DJ2", kThresholdDJ2}};
119  if(!fSelectedTriggers.size()) return false; // no trigger selected
120 
121  // Decode trigger clusters
122  const auto selclusters = GetAcceptedTriggerClusters(fInputEvent->GetFiredTriggerClasses().Data());
123 
124  auto findTriggerType = [](const std::vector<TString> &triggers, TString type) -> bool {
125  bool found = false;
126  for(const auto t : triggers) {
127  if(t.Contains(type)) {
128  found = true;
129  break;
130  }
131  }
132  return found;
133  };
134  bool isMB = std::find(fSelectedTriggers.begin(), fSelectedTriggers.end(), "MB") != fSelectedTriggers.end(),
135  isEG = findTriggerType(fSelectedTriggers, "EG"),
136  isDG = findTriggerType(fSelectedTriggers, "DG"),
137  isEJ = findTriggerType(fSelectedTriggers, "EJ"),
138  isDJ = findTriggerType(fSelectedTriggers, "DJ");
139 
140  std::vector<const AliEMCALTriggerPatchInfo *> EGApatches, DGApatches, EJEpatches, DJEpatches;
141  if(isMB || isEG) EGApatches = SelectAllPatchesByType(*fTriggerPatchInfo, kEGApatches);
142  if(isMB || isDG) DGApatches = SelectAllPatchesByType(*fTriggerPatchInfo, kDGApatches);
143  if(isMB || isEJ) EJEpatches = SelectAllPatchesByType(*fTriggerPatchInfo, kEJEpatches);
144  if(isMB || isDJ) DJEpatches = SelectAllPatchesByType(*fTriggerPatchInfo, kDJEpatches);
145 
146  for(const auto &t : fSelectedTriggers) {
147  if(!std::find(kNamesTriggerClasses.begin(), kNamesTriggerClasses.end(), t.Data())) continue;
148  if(t == "MB") {
149  // Min bias: Only fill patch ADC spectra all patches
150  for(auto patch : EGApatches) fHistos->FillTH1("hPatchADCEGAMB", patch->GetADCAmp());
151  for(auto patch : DGApatches) fHistos->FillTH1("hPatchADCDGAMB", patch->GetADCAmp());
152  for(auto patch : EJEpatches) fHistos->FillTH1("hPatchADCEJEMB", patch->GetADCAmp());
153  for(auto patch : DJEpatches) fHistos->FillTH1("hPatchADCDJEMB", patch->GetADCAmp());
154 
155  continue;
156  } else if(std::find(kNamesTriggerClasses.begin(), kNamesTriggerClasses.end(), t.Data()) != kNamesTriggerClasses.end()) {
157  const char detector = t[0];
158  const char *patchtype = ((t[1] == 'G') ? "GA" : "JE");
159  std::vector<const AliEMCALTriggerPatchInfo *> &patchhandler = (detector == 'E' ? (t[1] == 'G' ? EGApatches : EJEpatches) : (t[1] == 'G' ? DGApatches : DJEpatches));
160  auto firedpatches = SelectFiredPatchesByTrigger(*fTriggerPatchInfo, kPatchIndex.find(t.Data())->second);
161  auto patchareas = GetNumberNonOverlappingPatchAreas(firedpatches);
162  for(auto p : patchhandler){
163  double point[3] = {static_cast<double>(p->GetADCAmp()), static_cast<double>(firedpatches.size()), static_cast<double>(patchareas)};
164  for(const auto &kc : selclusters) {
165  fHistos->FillTH1(Form("hPatchADC%c%s%s%s", detector, patchtype, t.Data(), kc.data()), p->GetADCAmp());
166  fHistos->FillTHnSparse(Form("hAllPatches%c%s%s%s", detector, patchtype, t.Data(), kc.data()), point);
167  }
168  }
169  for(auto p : firedpatches) {
170  double point[5] = {static_cast<double>(p->GetADCAmp()), static_cast<double>(p->GetColStart()), static_cast<double>(p->GetRowStart()), static_cast<double>(firedpatches.size()), static_cast<double>(patchareas)};
171  for(const auto &kc : selclusters)
172  fHistos->FillTHnSparse(Form("hFiredPatches%c%s%s%s", detector, patchtype, t.Data(), kc.data()), point);
173  }
174  }
175  }
176  return true;
177 }
178 
179 std::vector<const AliEMCALTriggerPatchInfo *> AliAnalysisTaskEmcalRecalcPatchesRef::SelectAllPatchesByType(const TClonesArray &list, EPatchType_t patchtype) const {
180  std::vector<const AliEMCALTriggerPatchInfo *> result;
181  for(auto p : list){
182  AliEMCALTriggerPatchInfo *patch = static_cast<AliEMCALTriggerPatchInfo *>(p);
183  if(!patch->IsRecalc()) continue;
184  bool selected(true);
185  switch(patchtype){
186  case kEGApatches: if(patch->IsDCalPHOS() || !patch->IsGammaHighRecalc()) selected = false; break;
187  case kDGApatches: if(patch->IsEMCal() || !patch->IsGammaHighRecalc()) selected = false; break;
188  case kEJEpatches: if(patch->IsDCalPHOS() || !patch->IsJetHighRecalc()) selected = false; break;
189  case kDJEpatches: if(patch->IsEMCal() || !patch->IsJetHighRecalc()) selected = false; break;
190  };
191  if(selected) result.emplace_back(patch);
192  }
193  return result;
194 }
195 
196 std::vector<const AliEMCALTriggerPatchInfo *> AliAnalysisTaskEmcalRecalcPatchesRef::SelectFiredPatchesByTrigger(const TClonesArray &list, ETriggerThreshold_t trigger) const {
197  std::vector<const AliEMCALTriggerPatchInfo *> result;
198  EPatchType_t patchtype;
199  switch(trigger) {
200  case kThresholdEG1: patchtype = kEGApatches; break;
201  case kThresholdEG2: patchtype = kEGApatches; break;
202  case kThresholdDG1: patchtype = kDGApatches; break;
203  case kThresholdDG2: patchtype = kDGApatches; break;
204  case kThresholdEJ1: patchtype = kEJEpatches; break;
205  case kThresholdEJ2: patchtype = kEJEpatches; break;
206  case kThresholdDJ1: patchtype = kDJEpatches; break;
207  case kThresholdDJ2: patchtype = kDJEpatches; break;
208  default: return result; // unsupported patch type - return empty list
209  };
210  for(auto patch : SelectAllPatchesByType(list, patchtype)) {
211  if(patch->GetADCAmp() > fOnlineThresholds[trigger]) result.emplace_back(patch);
212  }
213  return result;
214 }
215 
216 std::vector<std::string> AliAnalysisTaskEmcalRecalcPatchesRef::GetAcceptedTriggerClusters(const char *triggerstring) const {
217  auto clusters = PWG::EMCAL::Triggerinfo::DecodeTriggerString(triggerstring);
218  std::vector<std::string> selclusters;
219  selclusters.push_back("ANY");
220  bool isCENT(false), isCENTNOTRD(false);
221  for(const auto &c : clusters){
222  if((c.Triggercluster() == "CENT") && !isCENT) { // don't count clusters multiple times
223  selclusters.push_back("CENT");
224  isCENT = true;
225  } else if((c.Triggercluster() == "CENTNOTRD") && !isCENTNOTRD) { // don't count clusters multiple times
226  selclusters.push_back("CENTNOTRD");
227  isCENTNOTRD = true;
228  }
229  }
230  if(isCENT && isCENTNOTRD) selclusters.push_back("BOTH");
231  else {
232  if(isCENT) selclusters.push_back("ONLYCENT");
233  if(isCENTNOTRD) selclusters.push_back("ONLYCENTNOTRD");
234  }
235  return selclusters;
236 }
237 
238 int AliAnalysisTaskEmcalRecalcPatchesRef::GetNumberNonOverlappingPatchAreas(const std::vector<const AliEMCALTriggerPatchInfo *> &firedpatches) const {
239  std::vector<const AliEMCALTriggerPatchInfo *> patchareas;
240  for(const auto patch : firedpatches) {
241  if(!patchareas.size()) patchareas.push_back(patch); // first patch always considered as new area
242  else {
243  bool overlapFound = false;
244  for(const auto refpatch : patchareas) {
245  if(!refpatch) {
246  AliErrorStream() << "Ref patch null" << std::endl;
247  AliErrorStream() << "Patchareas has size " << patchareas.size() << std::endl;
248  AliErrorStream() << "Firedpatches has size " << firedpatches.size() << std::endl;
249  }
250  if(!patch){
251  AliErrorStream() << "Test patch null" << std::endl;
252  AliErrorStream() << "Patchareas has size " << patchareas.size() << std::endl;
253  AliErrorStream() << "Firedpatches has size " << firedpatches.size() << std::endl;
254  }
255  if(HasOverlap(*refpatch, *patch)) {
256  // patch has overlap with allready accepted patch - discard
257  overlapFound = true;
258  break;
259  }
260  }
261  if(!overlapFound) patchareas.emplace_back(patch); // New non-overlapping patch found
262  }
263  }
264  return patchareas.size();
265 }
266 
267 bool AliAnalysisTaskEmcalRecalcPatchesRef::HasOverlap(const AliEMCALTriggerPatchInfo &ref, const AliEMCALTriggerPatchInfo &test) const {
268  int testcolmin = test.GetColStart(), testcolmax = test.GetColStart()+test.GetPatchSize()-1,
269  testrowmin = test.GetRowStart(), testrowmax = test.GetRowStart()+test.GetPatchSize()-1,
270  refcolmin = ref.GetColStart(), refcolmax = ref.GetColStart()+ref.GetPatchSize()-1,
271  refrowmin = ref.GetRowStart(), refrowmax = ref.GetRowStart()+ref.GetPatchSize()-1;
272  if((InRange(testcolmin, refcolmin, refcolmax) && InRange(testrowmin, refrowmin, refrowmax)) ||
273  (InRange(testcolmax, refcolmin, refcolmax) && InRange(testrowmax, refrowmin, refrowmax))) return true;
274  return false;
275 }
276 
278  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
279  if(!mgr){
280  std::cerr << "No analysis manager defined" << std::endl;
281  return nullptr;
282  }
283 
284  std::stringstream taskname;
285  taskname << "EmcalRecalcPatches_" << suffix;
287  mgr->AddTask(task);
288  task->SetEnableSumw2(true);
289 
298 
299  std::stringstream outfilename, outlistname;
300  outfilename << mgr->GetCommonFileName() << ":" << taskname.str();
301  outlistname << "Histos_" << taskname.str();
302  mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
303  mgr->ConnectOutput(task, 1, mgr->CreateContainer(outlistname.str().data(), TList::Class(), AliAnalysisManager::kOutputContainer, outfilename.str().data()));
304 
305  return task;
306 }
std::vector< TString > fSelectedTriggers
! Triggers selected for given event
bool HasOverlap(const AliEMCALTriggerPatchInfo &ref, const AliEMCALTriggerPatchInfo &test) const
Class creating a linear binning, used in the histogram manager.
TCanvas * c
Definition: TestFitELoss.C:172
Interface for binnings used by the histogram handler.
Definition: TBinning.h:21
static std::vector< PWG::EMCAL::Triggerinfo > DecodeTriggerString(const std::string &triggerstring)
Decoding trigger string.
void FillTHnSparse(const char *name, const double *x, double weight=1., Option_t *opt="")
static AliAnalysisTaskEmcalRecalcPatchesRef * AddTaskEmcalRecalcPatches(const char *suffix)
std::vector< const AliEMCALTriggerPatchInfo * > SelectFiredPatchesByTrigger(const TClonesArray &patches, ETriggerThreshold_t trigger) const
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
int GetNumberNonOverlappingPatchAreas(const std::vector< const AliEMCALTriggerPatchInfo * > &diredpatches) const
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
virtual bool Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
std::vector< std::string > GetAcceptedTriggerClusters(const char *triggerstring) const
Analysis of high- tracks in triggered events.
TClonesArray * fTriggerPatchInfo
!trigger patch info array
void test(int runnumber=195345)
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.
std::vector< const AliEMCALTriggerPatchInfo * > SelectAllPatchesByType(const TClonesArray &patches, EPatchType_t patchtype) const