AliPhysics  e6c8d43 (e6c8d43)
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 <TObjString.h>
38 #include <TString.h>
39 #include "AliAnalysisManager.h"
40 #include "AliEMCALTriggerPatchInfo.h"
43 
45 
46 using namespace EMCalTriggerPtAnalysis;
47 
50  fEnableSumw2(false),
51  fOnlineThresholds(),
52  fSwapPatches(false),
53  fRequiredOverlaps(),
54  fExcludedOverlaps()
55 {
56  SetCaloTriggerPatchInfoName("EmcalTriggers");
57 }
58 
61  fEnableSumw2(false),
63  fSwapPatches(false),
66 {
67  SetCaloTriggerPatchInfoName("EmcalTriggers");
68  fOnlineThresholds.Set(8);
69 }
70 
72  const std::array<std::string, 9> kNamesTriggerClasses = {{"MB", "EG1", "EG2", "DG1", "DG2", "EJ1", "EJ2", "DJ1", "DJ2"}};
73  const std::array<std::string, 6> kNamesTriggerClusters = {{"ANY", "CENT", "CENTNOTRD", "BOTH", "ONLYCENT", "ONLYCENTNOTRD"}};
74  const std::array<std::string, 4> kNamesPatchTypes = {{"EGA", "DGA", "EJE", "DJE"}};
75 
76  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);
77  const TBinning *firedpatchbinning[5] = {&ADCBinning, &colbinning, &rowbinning, &npatchbinning, &noverlapbinning},
78  *allpatchbinning[3] = {&ADCBinning, &npatchbinning, &noverlapbinning};
79  // Create event counters
80  for(const auto &kt : kNamesTriggerClasses){
81  // Min. Bias: Only CENT cluster - no distinction between trigger clusters
82  // EMCAL triggers: Distinction between CENT and CENTNOTRD necessary
83  if(kt.find("MB") != std::string::npos)
84  fHistos->CreateTH1(Form("hEventCounter%s", kt.data()), Form("Event counter for %s", kt.data()), 1, 0.5, 1.5);
85  else {
86  for(const auto &kc : kNamesTriggerClusters) {
87  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);
88  }
89  }
90  }
91 
92  // Min. Bias: Create patch energy spectra (all patches) for EMCAL and DCAL
93  // Min. Bias trigger only in CENT cluster
94  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" : "");
95 
96  // Triggers: Create trigger spectra and THnSparse of firing patches
97  for(const auto &kt : kNamesTriggerClasses) {
98  if(kt == "MB") continue;
99  const char detector = kt[0];
100  const char *patchtype = ((kt[1] == 'G') ? "GA" : "JE");
101  // distinction between trigger clusters
102  for(const auto &kc : kNamesTriggerClusters){
103  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" : "");
104  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" : "");
105  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" : "");
106  }
107  }
108 }
109 
111  // handle overlaps
112  if(fRequiredOverlaps.GetEntries() || fExcludedOverlaps.GetEntries()){
113  if(fRequiredOverlaps.GetEntries()){
114  Bool_t allFound(true);
115  for(auto t : fRequiredOverlaps){
116  auto trgstr = static_cast<TObjString *>(t)->String();
117  if(std::find_if(fSelectedTriggers.begin(), fSelectedTriggers.end(), [&trgstr](const TString &seltrigger) -> bool { return seltrigger.Contains(trgstr); }) == fSelectedTriggers.end()) {
118  allFound = false;
119  break;
120  }
121  }
122  if(!allFound) return false;
123  }
124 
125  if(fExcludedOverlaps.GetEntries()){
126  Bool_t oneFound(false);
127  for(auto t : fExcludedOverlaps){
128  auto trgstr = static_cast<TObjString *>(t)->String();
129  if(std::find_if(fSelectedTriggers.begin(), fSelectedTriggers.end(), [&trgstr](const TString &seltrigger) -> bool { return seltrigger.Contains(trgstr); }) != fSelectedTriggers.end()) {
130  oneFound = true;
131  break;
132  }
133  }
134  if(oneFound) return false;
135  }
136  }
137  return true;
138 }
139 
141  const std::array<std::string, 9> kNamesTriggerClasses = {{"MB", "EG1", "EG2", "DG1", "DG2", "EJ1", "EJ2", "DJ1", "DJ2"}};
142  const auto selclusters = GetAcceptedTriggerClusters(fInputEvent->GetFiredTriggerClasses().Data());
143 
144  auto findTriggerType = [](const std::vector<TString> &triggers, TString type) -> bool {
145  bool found = false;
146  for(const auto t : triggers) {
147  if(t.Contains(type)) {
148  found = true;
149  break;
150  }
151  }
152  return found;
153  };
154 
155  std::vector<TString> handledtriggers;
156  for(auto t : kNamesTriggerClasses) {
157  if(findTriggerType(fSelectedTriggers, t.data())){
158  handledtriggers.emplace_back(t);
159  }
160  }
161 
162  for(const auto &kt : handledtriggers) {
163  if(kt == "MB") fHistos->FillTH1(Form("hEventCounter%s", kt.Data()), 1.);
164  else {
165  for(const auto &kc : selclusters) fHistos->FillTH1(Form("hEventCounter%s%s", kt.Data(), kc.data()), 1.);
166  }
167  }
168 }
169 
171  const std::array<std::string, 9> kNamesTriggerClasses = {{"MB", "EG1", "EG2", "DG1", "DG2", "EJ1", "EJ2", "DJ1", "DJ2"}};
172  const std::unordered_map<std::string, ETriggerThreshold_t> kPatchIndex = {{"EG1", kThresholdEG1},{"EG2", kThresholdEG2},{"DG1", kThresholdDG1},{"DG2", kThresholdDG2},
173  {"EJ1", kThresholdEJ1},{"EJ2", kThresholdEJ2},{"DJ1", kThresholdDJ1},{"DJ2", kThresholdDJ2}};
174  if(!fSelectedTriggers.size()) return false; // no trigger selected
175  AliDebugStream(1) << "Found triggers" << std::endl;
176  for(auto t : fSelectedTriggers) AliDebugStream(1) << t << std::endl;
177  AliDebugStream(1) << "Trigger patch container has " << fTriggerPatchInfo->GetEntries() << " patches" << std::endl;
178 
179  // Decode trigger clusters
180  const auto selclusters = GetAcceptedTriggerClusters(fInputEvent->GetFiredTriggerClasses().Data());
181 
182  auto findTriggerType = [](const std::vector<TString> &triggers, TString type) -> bool {
183  bool found = false;
184  for(const auto t : triggers) {
185  if(t.Contains(type)) {
186  found = true;
187  break;
188  }
189  }
190  return found;
191  };
192  bool isMB = std::find(fSelectedTriggers.begin(), fSelectedTriggers.end(), "MB") != fSelectedTriggers.end(),
193  isEG = findTriggerType(fSelectedTriggers, "EG"),
194  isDG = findTriggerType(fSelectedTriggers, "DG"),
195  isEJ = findTriggerType(fSelectedTriggers, "EJ"),
196  isDJ = findTriggerType(fSelectedTriggers, "DJ");
197 
198  std::vector<TString> handledtriggers;
199  for(auto t : kNamesTriggerClasses) {
200  if(findTriggerType(fSelectedTriggers, t.data())){
201  handledtriggers.emplace_back(t);
202  }
203  }
204 
205  AliDebugStream(1) << "Processing triggers" << std::endl;
206  for(auto t : handledtriggers) AliDebugStream(1) << t << std::endl;
207  if(!handledtriggers.size()){
208  AliDebugStream(1) << "No supported trigger class found " << std::endl;
209  return false;
210  }
211 
212  std::vector<const AliEMCALTriggerPatchInfo *> EGApatches, DGApatches, EJEpatches, DJEpatches;
213  if(isMB || isEG) EGApatches = SelectAllPatchesByType(*fTriggerPatchInfo, fSwapPatches ? kEJEpatches : kEGApatches);
214  if(isMB || isDG) DGApatches = SelectAllPatchesByType(*fTriggerPatchInfo, fSwapPatches ? kDJEpatches : kDGApatches);
215  if(isMB || isEJ) EJEpatches = SelectAllPatchesByType(*fTriggerPatchInfo, fSwapPatches ? kEGApatches : kEJEpatches);
216  if(isMB || isDJ) DJEpatches = SelectAllPatchesByType(*fTriggerPatchInfo, fSwapPatches ? kDGApatches : kDJEpatches);
217 
218  for(const auto &t : handledtriggers) {
219  if(t == "MB") {
220  // Min bias: Only fill patch ADC spectra all patches
221  for(auto patch : EGApatches) fHistos->FillTH1("hPatchADCEGAMB", patch->GetADCAmp());
222  for(auto patch : DGApatches) fHistos->FillTH1("hPatchADCDGAMB", patch->GetADCAmp());
223  for(auto patch : EJEpatches) fHistos->FillTH1("hPatchADCEJEMB", patch->GetADCAmp());
224  for(auto patch : DJEpatches) fHistos->FillTH1("hPatchADCDJEMB", patch->GetADCAmp());
225  } else {
226  const char detector = t[0];
227  const char *patchtype = ((t[1] == 'G') ? "GA" : "JE");
228  std::vector<const AliEMCALTriggerPatchInfo *> &patchhandler = (detector == 'E' ? (t[1] == 'G' ? EGApatches : EJEpatches) : (t[1] == 'G' ? DGApatches : DJEpatches));
229  auto firedpatches = SelectFiredPatchesByTrigger(*fTriggerPatchInfo, kPatchIndex.find(t.Data())->second);
230  auto patchareas = GetNumberNonOverlappingPatchAreas(firedpatches);
231  AliDebugStream(3) << "Trigger " << t << ", patches " << patchhandler.size() << ", firing " << firedpatches.size() << std::endl;
232  for(auto p : patchhandler){
233  double point[3] = {static_cast<double>(p->GetADCAmp()), static_cast<double>(firedpatches.size()), static_cast<double>(patchareas)};
234  for(const auto &kc : selclusters) {
235  fHistos->FillTH1(Form("hPatchADC%c%s%s%s", detector, patchtype, t.Data(), kc.data()), p->GetADCAmp());
236  fHistos->FillTHnSparse(Form("hAllPatches%c%s%s%s", detector, patchtype, t.Data(), kc.data()), point);
237  }
238  }
239  for(auto p : firedpatches) {
240  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)};
241  for(const auto &kc : selclusters)
242  fHistos->FillTHnSparse(Form("hFiredPatches%c%s%s%s", detector, patchtype, t.Data(), kc.data()), point);
243  }
244  }
245  }
246  return true;
247 }
248 
249 std::vector<const AliEMCALTriggerPatchInfo *> AliAnalysisTaskEmcalRecalcPatchesRef::SelectAllPatchesByType(const TClonesArray &list, EPatchType_t patchtype) const {
250  AliDebugStream(2) << "Selecting all patches for trigger " << static_cast<int>(patchtype) << std::endl;
251  std::vector<const AliEMCALTriggerPatchInfo *> result;
252  for(auto p : list){
253  AliEMCALTriggerPatchInfo *patch = static_cast<AliEMCALTriggerPatchInfo *>(p);
254  if(!patch->IsRecalc()) continue;
255  bool selected(true);
256  switch(patchtype){
257  case kEGApatches: if(patch->IsDCalPHOS() || !patch->IsGammaHighRecalc()) selected = false; break;
258  case kDGApatches: if(patch->IsEMCal() || !patch->IsGammaHighRecalc()) selected = false; break;
259  case kEJEpatches: if(patch->IsDCalPHOS() || !patch->IsJetHighRecalc()) selected = false; break;
260  case kDJEpatches: if(patch->IsEMCal() || !patch->IsJetHighRecalc()) selected = false; break;
261  };
262  if(selected) result.emplace_back(patch);
263  }
264  AliDebugStream(2) << "In: " << list.GetEntries() << ", out: " << result.size() << std::endl;
265  return result;
266 }
267 
268 std::vector<const AliEMCALTriggerPatchInfo *> AliAnalysisTaskEmcalRecalcPatchesRef::SelectFiredPatchesByTrigger(const TClonesArray &list, ETriggerThreshold_t trigger) const {
269  std::vector<const AliEMCALTriggerPatchInfo *> result;
270  EPatchType_t patchtype;
271  switch(trigger) {
272  case kThresholdEG1: patchtype = kEGApatches; break;
273  case kThresholdEG2: patchtype = kEGApatches; break;
274  case kThresholdDG1: patchtype = kDGApatches; break;
275  case kThresholdDG2: patchtype = kDGApatches; break;
276  case kThresholdEJ1: patchtype = kEJEpatches; break;
277  case kThresholdEJ2: patchtype = kEJEpatches; break;
278  case kThresholdDJ1: patchtype = kDJEpatches; break;
279  case kThresholdDJ2: patchtype = kDJEpatches; break;
280  default: return result; // unsupported patch type - return empty list
281  };
282  for(auto patch : SelectAllPatchesByType(list, patchtype)) {
283  if(patch->GetADCAmp() > fOnlineThresholds[trigger]) result.emplace_back(patch);
284  }
285  return result;
286 }
287 
288 std::vector<std::string> AliAnalysisTaskEmcalRecalcPatchesRef::GetAcceptedTriggerClusters(const char *triggerstring) const {
289  auto clusters = PWG::EMCAL::Triggerinfo::DecodeTriggerString(triggerstring);
290  std::vector<std::string> selclusters;
291  selclusters.push_back("ANY");
292  bool isCENT(false), isCENTNOTRD(false);
293  for(const auto &c : clusters){
294  if((c.Triggercluster() == "CENT") && !isCENT) { // don't count clusters multiple times
295  selclusters.push_back("CENT");
296  isCENT = true;
297  } else if((c.Triggercluster() == "CENTNOTRD") && !isCENTNOTRD) { // don't count clusters multiple times
298  selclusters.push_back("CENTNOTRD");
299  isCENTNOTRD = true;
300  }
301  }
302  if(isCENT && isCENTNOTRD) selclusters.push_back("BOTH");
303  else {
304  if(isCENT) selclusters.push_back("ONLYCENT");
305  if(isCENTNOTRD) selclusters.push_back("ONLYCENTNOTRD");
306  }
307  return selclusters;
308 }
309 
310 int AliAnalysisTaskEmcalRecalcPatchesRef::GetNumberNonOverlappingPatchAreas(const std::vector<const AliEMCALTriggerPatchInfo *> &firedpatches) const {
311  std::vector<const AliEMCALTriggerPatchInfo *> patchareas;
312  for(const auto patch : firedpatches) {
313  if(!patchareas.size()) patchareas.push_back(patch); // first patch always considered as new area
314  else {
315  bool overlapFound = false;
316  for(const auto refpatch : patchareas) {
317  if(!refpatch) {
318  AliErrorStream() << "Ref patch null" << std::endl;
319  AliErrorStream() << "Patchareas has size " << patchareas.size() << std::endl;
320  AliErrorStream() << "Firedpatches has size " << firedpatches.size() << std::endl;
321  }
322  if(!patch){
323  AliErrorStream() << "Test patch null" << std::endl;
324  AliErrorStream() << "Patchareas has size " << patchareas.size() << std::endl;
325  AliErrorStream() << "Firedpatches has size " << firedpatches.size() << std::endl;
326  }
327  if(HasOverlap(*refpatch, *patch)) {
328  // patch has overlap with allready accepted patch - discard
329  overlapFound = true;
330  break;
331  }
332  }
333  if(!overlapFound) patchareas.emplace_back(patch); // New non-overlapping patch found
334  }
335  }
336  return patchareas.size();
337 }
338 
339 bool AliAnalysisTaskEmcalRecalcPatchesRef::HasOverlap(const AliEMCALTriggerPatchInfo &ref, const AliEMCALTriggerPatchInfo &test) const {
340  int testcolmin = test.GetColStart(), testcolmax = test.GetColStart()+test.GetPatchSize()-1,
341  testrowmin = test.GetRowStart(), testrowmax = test.GetRowStart()+test.GetPatchSize()-1,
342  refcolmin = ref.GetColStart(), refcolmax = ref.GetColStart()+ref.GetPatchSize()-1,
343  refrowmin = ref.GetRowStart(), refrowmax = ref.GetRowStart()+ref.GetPatchSize()-1;
344  if((InRange(testcolmin, refcolmin, refcolmax) && InRange(testrowmin, refrowmin, refrowmax)) ||
345  (InRange(testcolmax, refcolmin, refcolmax) && InRange(testrowmax, refrowmin, refrowmax))) return true;
346  return false;
347 }
348 
350  if(fRequiredOverlaps.FindObject(trigger)) return;
351  fRequiredOverlaps.Add(new TObjString(trigger));
352 }
353 
355  if(fExcludedOverlaps.FindObject(trigger)) return;
356  fExcludedOverlaps.Add(new TObjString(trigger));
357 }
358 
360  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
361  if(!mgr){
362  std::cerr << "No analysis manager defined" << std::endl;
363  return nullptr;
364  }
365 
366  std::stringstream taskname;
367  taskname << "EmcalRecalcPatches_" << suffix;
369  mgr->AddTask(task);
370  task->SetEnableSumw2(true);
371 
380 
381  std::stringstream outfilename, outlistname;
382  outfilename << mgr->GetCommonFileName() << ":" << taskname.str();
383  outlistname << "Histos_" << taskname.str();
384  mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
385  mgr->ConnectOutput(task, 1, mgr->CreateContainer(outlistname.str().data(), TList::Class(), AliAnalysisManager::kOutputContainer, outfilename.str().data()));
386 
387  return task;
388 }
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 SetCaloTriggerPatchInfoName(const char *n)
TObjArray fRequiredOverlaps
Add option to require overlap with certain triggers.
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.
TObjArray fExcludedOverlaps
Add option to exclude overlap with certain triggers.
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)
bool Bool_t
Definition: External.C:53
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