AliPhysics  ced2227 (ced2227)
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 
106  auto findTriggerType = [](const std::vector<TString> &triggers, TString type) -> bool {
107  bool found = false;
108  for(const auto t : triggers) {
109  if(t.Contains(type)) {
110  found = true;
111  break;
112  }
113  }
114  return found;
115  };
116 
117  std::vector<TString> handledtriggers;
118  for(auto t : kNamesTriggerClasses) {
119  if(findTriggerType(fSelectedTriggers, t.data())){
120  handledtriggers.emplace_back(t);
121  }
122  }
123 
124  for(const auto &kt : handledtriggers) {
125  if(kt == "MB") fHistos->FillTH1(Form("hEventCounter%s", kt.Data()), 1.);
126  else {
127  for(const auto &kc : selclusters) fHistos->FillTH1(Form("hEventCounter%s%s", kt.Data(), kc.data()), 1.);
128  }
129  }
130 }
131 
133  const std::array<std::string, 9> kNamesTriggerClasses = {{"MB", "EG1", "EG2", "DG1", "DG2", "EJ1", "EJ2", "DJ1", "DJ2"}};
134  const std::unordered_map<std::string, ETriggerThreshold_t> kPatchIndex = {{"EG1", kThresholdEG1},{"EG2", kThresholdEG2},{"DG1", kThresholdDG1},{"DG2", kThresholdDG2},
135  {"EJ1", kThresholdEJ1},{"EJ2", kThresholdEJ2},{"DJ1", kThresholdDJ1},{"DJ2", kThresholdDJ2}};
136  if(!fSelectedTriggers.size()) return false; // no trigger selected
137  AliDebugStream(1) << "Found triggers" << std::endl;
138  for(auto t : fSelectedTriggers) AliDebugStream(1) << t << std::endl;
139  AliDebugStream(1) << "Trigger patch container has " << fTriggerPatchInfo->GetEntries() << " patches" << std::endl;
140 
141  // Decode trigger clusters
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  bool isMB = std::find(fSelectedTriggers.begin(), fSelectedTriggers.end(), "MB") != fSelectedTriggers.end(),
155  isEG = findTriggerType(fSelectedTriggers, "EG"),
156  isDG = findTriggerType(fSelectedTriggers, "DG"),
157  isEJ = findTriggerType(fSelectedTriggers, "EJ"),
158  isDJ = findTriggerType(fSelectedTriggers, "DJ");
159 
160  std::vector<TString> handledtriggers;
161  for(auto t : kNamesTriggerClasses) {
162  if(findTriggerType(fSelectedTriggers, t.data())){
163  handledtriggers.emplace_back(t);
164  }
165  }
166 
167  AliDebugStream(1) << "Processing triggers" << std::endl;
168  for(auto t : handledtriggers) AliDebugStream(1) << t << std::endl;
169  if(!handledtriggers.size()){
170  AliDebugStream(1) << "No supported trigger class found " << std::endl;
171  return false;
172  }
173 
174  std::vector<const AliEMCALTriggerPatchInfo *> EGApatches, DGApatches, EJEpatches, DJEpatches;
175  if(isMB || isEG) EGApatches = SelectAllPatchesByType(*fTriggerPatchInfo, kEGApatches);
176  if(isMB || isDG) DGApatches = SelectAllPatchesByType(*fTriggerPatchInfo, kDGApatches);
177  if(isMB || isEJ) EJEpatches = SelectAllPatchesByType(*fTriggerPatchInfo, kEJEpatches);
178  if(isMB || isDJ) DJEpatches = SelectAllPatchesByType(*fTriggerPatchInfo, kDJEpatches);
179 
180  for(const auto &t : handledtriggers) {
181  if(t == "MB") {
182  // Min bias: Only fill patch ADC spectra all patches
183  for(auto patch : EGApatches) fHistos->FillTH1("hPatchADCEGAMB", patch->GetADCAmp());
184  for(auto patch : DGApatches) fHistos->FillTH1("hPatchADCDGAMB", patch->GetADCAmp());
185  for(auto patch : EJEpatches) fHistos->FillTH1("hPatchADCEJEMB", patch->GetADCAmp());
186  for(auto patch : DJEpatches) fHistos->FillTH1("hPatchADCDJEMB", patch->GetADCAmp());
187  } else {
188  const char detector = t[0];
189  const char *patchtype = ((t[1] == 'G') ? "GA" : "JE");
190  std::vector<const AliEMCALTriggerPatchInfo *> &patchhandler = (detector == 'E' ? (t[1] == 'G' ? EGApatches : EJEpatches) : (t[1] == 'G' ? DGApatches : DJEpatches));
191  auto firedpatches = SelectFiredPatchesByTrigger(*fTriggerPatchInfo, kPatchIndex.find(t.Data())->second);
192  auto patchareas = GetNumberNonOverlappingPatchAreas(firedpatches);
193  AliDebugStream(3) << "Trigger " << t << ", patches " << patchhandler.size() << ", firing " << firedpatches.size() << std::endl;
194  for(auto p : patchhandler){
195  double point[3] = {static_cast<double>(p->GetADCAmp()), static_cast<double>(firedpatches.size()), static_cast<double>(patchareas)};
196  for(const auto &kc : selclusters) {
197  fHistos->FillTH1(Form("hPatchADC%c%s%s%s", detector, patchtype, t.Data(), kc.data()), p->GetADCAmp());
198  fHistos->FillTHnSparse(Form("hAllPatches%c%s%s%s", detector, patchtype, t.Data(), kc.data()), point);
199  }
200  }
201  for(auto p : firedpatches) {
202  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)};
203  for(const auto &kc : selclusters)
204  fHistos->FillTHnSparse(Form("hFiredPatches%c%s%s%s", detector, patchtype, t.Data(), kc.data()), point);
205  }
206  }
207  }
208  return true;
209 }
210 
211 std::vector<const AliEMCALTriggerPatchInfo *> AliAnalysisTaskEmcalRecalcPatchesRef::SelectAllPatchesByType(const TClonesArray &list, EPatchType_t patchtype) const {
212  AliDebugStream(2) << "Selecting all patches for trigger " << static_cast<int>(patchtype) << std::endl;
213  std::vector<const AliEMCALTriggerPatchInfo *> result;
214  for(auto p : list){
215  AliEMCALTriggerPatchInfo *patch = static_cast<AliEMCALTriggerPatchInfo *>(p);
216  if(!patch->IsRecalc()) continue;
217  bool selected(true);
218  switch(patchtype){
219  case kEGApatches: if(patch->IsDCalPHOS() || !patch->IsGammaHighRecalc()) selected = false; break;
220  case kDGApatches: if(patch->IsEMCal() || !patch->IsGammaHighRecalc()) selected = false; break;
221  case kEJEpatches: if(patch->IsDCalPHOS() || !patch->IsJetHighRecalc()) selected = false; break;
222  case kDJEpatches: if(patch->IsEMCal() || !patch->IsJetHighRecalc()) selected = false; break;
223  };
224  if(selected) result.emplace_back(patch);
225  }
226  AliDebugStream(2) << "In: " << list.GetEntries() << ", out: " << result.size() << std::endl;
227  return result;
228 }
229 
230 std::vector<const AliEMCALTriggerPatchInfo *> AliAnalysisTaskEmcalRecalcPatchesRef::SelectFiredPatchesByTrigger(const TClonesArray &list, ETriggerThreshold_t trigger) const {
231  std::vector<const AliEMCALTriggerPatchInfo *> result;
232  EPatchType_t patchtype;
233  switch(trigger) {
234  case kThresholdEG1: patchtype = kEGApatches; break;
235  case kThresholdEG2: patchtype = kEGApatches; break;
236  case kThresholdDG1: patchtype = kDGApatches; break;
237  case kThresholdDG2: patchtype = kDGApatches; break;
238  case kThresholdEJ1: patchtype = kEJEpatches; break;
239  case kThresholdEJ2: patchtype = kEJEpatches; break;
240  case kThresholdDJ1: patchtype = kDJEpatches; break;
241  case kThresholdDJ2: patchtype = kDJEpatches; break;
242  default: return result; // unsupported patch type - return empty list
243  };
244  for(auto patch : SelectAllPatchesByType(list, patchtype)) {
245  if(patch->GetADCAmp() > fOnlineThresholds[trigger]) result.emplace_back(patch);
246  }
247  return result;
248 }
249 
250 std::vector<std::string> AliAnalysisTaskEmcalRecalcPatchesRef::GetAcceptedTriggerClusters(const char *triggerstring) const {
251  auto clusters = PWG::EMCAL::Triggerinfo::DecodeTriggerString(triggerstring);
252  std::vector<std::string> selclusters;
253  selclusters.push_back("ANY");
254  bool isCENT(false), isCENTNOTRD(false);
255  for(const auto &c : clusters){
256  if((c.Triggercluster() == "CENT") && !isCENT) { // don't count clusters multiple times
257  selclusters.push_back("CENT");
258  isCENT = true;
259  } else if((c.Triggercluster() == "CENTNOTRD") && !isCENTNOTRD) { // don't count clusters multiple times
260  selclusters.push_back("CENTNOTRD");
261  isCENTNOTRD = true;
262  }
263  }
264  if(isCENT && isCENTNOTRD) selclusters.push_back("BOTH");
265  else {
266  if(isCENT) selclusters.push_back("ONLYCENT");
267  if(isCENTNOTRD) selclusters.push_back("ONLYCENTNOTRD");
268  }
269  return selclusters;
270 }
271 
272 int AliAnalysisTaskEmcalRecalcPatchesRef::GetNumberNonOverlappingPatchAreas(const std::vector<const AliEMCALTriggerPatchInfo *> &firedpatches) const {
273  std::vector<const AliEMCALTriggerPatchInfo *> patchareas;
274  for(const auto patch : firedpatches) {
275  if(!patchareas.size()) patchareas.push_back(patch); // first patch always considered as new area
276  else {
277  bool overlapFound = false;
278  for(const auto refpatch : patchareas) {
279  if(!refpatch) {
280  AliErrorStream() << "Ref patch null" << std::endl;
281  AliErrorStream() << "Patchareas has size " << patchareas.size() << std::endl;
282  AliErrorStream() << "Firedpatches has size " << firedpatches.size() << std::endl;
283  }
284  if(!patch){
285  AliErrorStream() << "Test patch null" << std::endl;
286  AliErrorStream() << "Patchareas has size " << patchareas.size() << std::endl;
287  AliErrorStream() << "Firedpatches has size " << firedpatches.size() << std::endl;
288  }
289  if(HasOverlap(*refpatch, *patch)) {
290  // patch has overlap with allready accepted patch - discard
291  overlapFound = true;
292  break;
293  }
294  }
295  if(!overlapFound) patchareas.emplace_back(patch); // New non-overlapping patch found
296  }
297  }
298  return patchareas.size();
299 }
300 
301 bool AliAnalysisTaskEmcalRecalcPatchesRef::HasOverlap(const AliEMCALTriggerPatchInfo &ref, const AliEMCALTriggerPatchInfo &test) const {
302  int testcolmin = test.GetColStart(), testcolmax = test.GetColStart()+test.GetPatchSize()-1,
303  testrowmin = test.GetRowStart(), testrowmax = test.GetRowStart()+test.GetPatchSize()-1,
304  refcolmin = ref.GetColStart(), refcolmax = ref.GetColStart()+ref.GetPatchSize()-1,
305  refrowmin = ref.GetRowStart(), refrowmax = ref.GetRowStart()+ref.GetPatchSize()-1;
306  if((InRange(testcolmin, refcolmin, refcolmax) && InRange(testrowmin, refrowmin, refrowmax)) ||
307  (InRange(testcolmax, refcolmin, refcolmax) && InRange(testrowmax, refrowmin, refrowmax))) return true;
308  return false;
309 }
310 
312  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
313  if(!mgr){
314  std::cerr << "No analysis manager defined" << std::endl;
315  return nullptr;
316  }
317 
318  std::stringstream taskname;
319  taskname << "EmcalRecalcPatches_" << suffix;
321  mgr->AddTask(task);
322  task->SetEnableSumw2(true);
323 
332 
333  std::stringstream outfilename, outlistname;
334  outfilename << mgr->GetCommonFileName() << ":" << taskname.str();
335  outlistname << "Histos_" << taskname.str();
336  mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
337  mgr->ConnectOutput(task, 1, mgr->CreateContainer(outlistname.str().data(), TList::Class(), AliAnalysisManager::kOutputContainer, outfilename.str().data()));
338 
339  return task;
340 }
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