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