AliPhysics  ced2227 (ced2227)
AliAnalysisTaskEmcalClusterNoiseTriggers.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 "AliAODInputHandler.h"
40 #include "AliClusterContainer.h"
41 #include "AliEMCALTriggerPatchInfo.h"
45 
47 
48 using namespace EMCalTriggerPtAnalysis;
49 
52  fEnableSumw2(false),
53  fOnlineThresholds()
54 {
55 
56 }
57 
60  fEnableSumw2(false),
62 {
63  fOnlineThresholds.Set(8);
64 }
65 
67  const std::array<std::string, 9> kNamesTriggerClasses = {{"MB", "EG1", "EG2", "DG1", "DG2", "EJ1", "EJ2", "DJ1", "DJ2"}};
68  const std::array<std::string, 6> kNamesTriggerClusters = {{"ANY", "CENT", "CENTNOTRD", "BOTH", "ONLYCENT", "ONLYCENTNOTRD"}};
69  const std::array<std::string, 2> kNamesDetectors = {{"EMCAL", "DCAL"}};
70 
71  TLinearBinning energybinning(2000, 0., 2000.), npatchbinning(51, -0.5, 50.5), noverlapbinning(21, -0.5, 20.5);
72  const TBinning *corrbinning[3] = {&energybinning, &npatchbinning, &noverlapbinning};
73  // Create event counters
74  for(const auto &kt : kNamesTriggerClasses){
75  // Min. Bias: Only CENT cluster - no distinction between trigger clusters
76  // EMCAL triggers: Distinction between CENT and CENTNOTRD necessary
77  if(kt.find("MB") != std::string::npos)
78  fHistos->CreateTH1(Form("hEventCounter%s", kt.data()), Form("Event counter for %s", kt.data()), 1, 0.5, 1.5);
79  else {
80  for(const auto &kc : kNamesTriggerClusters) {
81  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);
82  }
83  }
84  }
85 
86  // Min. Bias: Create patch energy spectra (all patches) for EMCAL and DCAL
87  // Min. Bias trigger only in CENT cluster
88  for(const auto &kd : kNamesDetectors) fHistos->CreateTH1(Form("hClusterEnergy%sMB", kd.data()), Form("%s cluster energy spectrum in MB events", kd.data()), 200, 0., 200., fEnableSumw2 ? "s" : "");
89 
90  // Triggers: Create trigger spectra and THnSparse of firing patches
91  for(const auto &kt : kNamesTriggerClasses) {
92  if(kt == "MB") continue;
93  // distinction between trigger clusters
94  for(const auto &kc : kNamesTriggerClusters){
95  fHistos->CreateTH1(Form("hClusterEnergy%s%s", kt.data(), kc.data()), Form("Cluster energy spectra in %s events (cluster %s)", kt.data(), kc.data()), 200, 0., 200., fEnableSumw2 ? "s" : "");
96  fHistos->CreateTHnSparse(Form("hClusterEnergyCorr%s%s", kt.data(), kc.data()), Form("Cluster energy spectra for trigger %s (cluster %s)", kt.data(), kc.data()), 3, corrbinning, fEnableSumw2 ? "s" : "");
97  }
98  }
99 }
100 
102  const std::array<std::string, 9> kNamesTriggerClasses = {{"MB", "EG1", "EG2", "DG1", "DG2", "EJ1", "EJ2", "DJ1", "DJ2"}};
103  const auto selclusters = GetAcceptedTriggerClusters(fInputEvent->GetFiredTriggerClasses().Data());
104  for(const auto &kt : kNamesTriggerClasses) {
105  if(std::find(fSelectedTriggers.begin(), fSelectedTriggers.end(), kt) != fSelectedTriggers.end()) {
106  if(kt == "MB") fHistos->FillTH1(Form("hEventCounter%s", kt.data()), 1.);
107  else {
108  for(const auto &kc : selclusters) fHistos->FillTH1(Form("hEventCounter%s%s", kt.data(), kc.data()), 1.);
109  }
110  }
111  }
112 }
113 
115  const std::array<std::string, 9> kNamesTriggerClasses = {{"MB", "EG1", "EG2", "DG1", "DG2", "EJ1", "EJ2", "DJ1", "DJ2"}};
116  const std::unordered_map<std::string, ETriggerThreshold_t> kPatchIndex = {{"EG1", kThresholdEG1},{"EG2", kThresholdEG2},{"DG1", kThresholdDG1},{"DG2", kThresholdDG2},
117  {"EJ1", kThresholdEJ1},{"EJ2", kThresholdEJ2},{"DJ1", kThresholdDJ1},{"DJ2", kThresholdDJ2}};
118  if(!fSelectedTriggers.size()) return false; // no trigger selected
119 
120  // Decode trigger clusters
121  const auto selclusters = GetAcceptedTriggerClusters(fInputEvent->GetFiredTriggerClasses().Data());
122 
123  auto findTriggerType = [](const std::vector<TString> &triggers, TString type) -> bool {
124  bool found = false;
125  for(const auto t : triggers) {
126  if(t.Contains(type)) {
127  found = true;
128  break;
129  }
130  }
131  return found;
132  };
133  bool isMB = std::find(fSelectedTriggers.begin(), fSelectedTriggers.end(), "MB") != fSelectedTriggers.end(),
134  isEG = findTriggerType(fSelectedTriggers, "EG"),
135  isDG = findTriggerType(fSelectedTriggers, "DG"),
136  isEJ = findTriggerType(fSelectedTriggers, "EJ"),
137  isDJ = findTriggerType(fSelectedTriggers, "DJ");
138 
139  auto clustercont = this->GetClusterContainer(0);
140 
141  for(const auto &t : fSelectedTriggers) {
142  if(!std::find(kNamesTriggerClasses.begin(), kNamesTriggerClasses.end(), t.Data())) continue;
143  if(t == "MB") {
144  // Min bias: Only fill patch ADC spectra all patches
145  for(auto c : clustercont->accepted()){
146  bool emccluster = (c->GetCellAbsId(0)<12288);
147  if(emccluster) fHistos->FillTH1("hClusterEnergyEMCALMB", c->GetNonLinCorrEnergy());
148  else fHistos->FillTH1("hClusterEnergyDCALMB", c->GetNonLinCorrEnergy());
149  }
150  } else if(std::find(kNamesTriggerClasses.begin(), kNamesTriggerClasses.end(), t.Data()) != kNamesTriggerClasses.end()) {
151  const char detector = t[0];
152  auto firedpatches = SelectFiredPatchesByTrigger(*fTriggerPatchInfo, kPatchIndex.find(t.Data())->second);
153  auto patchareas = GetNumberNonOverlappingPatchAreas(firedpatches);
154  for(auto c : clustercont->accepted()){
155  bool emccluster = (c->GetCellAbsId(0)<12288);
156  if((detector == 'E' && !emccluster) || (detector == 'D' && emccluster)) continue;
157  double point[3] = {c->GetNonLinCorrEnergy(), static_cast<double>(firedpatches.size()), static_cast<double>(patchareas)};
158  for(const auto &kc : selclusters) {
159  fHistos->FillTH1(Form("hClusterEnergy%s%s", t.Data(), kc.data()), c->GetNonLinCorrEnergy());
160  fHistos->FillTHnSparse(Form("hClusterEnergyCorr%s%s", t.Data(), kc.data()), point);
161  }
162  }
163  }
164  }
165  return true;
166 }
167 
168 std::vector<const AliEMCALTriggerPatchInfo *> AliAnalysisTaskEmcalClusterNoiseTriggers::SelectAllPatchesByType(const TClonesArray &list, EPatchType_t patchtype) const {
169  std::vector<const AliEMCALTriggerPatchInfo *> result;
170  for(auto p : list){
171  AliEMCALTriggerPatchInfo *patch = static_cast<AliEMCALTriggerPatchInfo *>(p);
172  if(!patch->IsRecalc()) continue;
173  bool selected(true);
174  switch(patchtype){
175  case kEGApatches: if(patch->IsDCalPHOS() || !patch->IsGammaHighRecalc()) selected = false; break;
176  case kDGApatches: if(patch->IsEMCal() || !patch->IsGammaHighRecalc()) selected = false; break;
177  case kEJEpatches: if(patch->IsDCalPHOS() || !patch->IsJetHighRecalc()) selected = false; break;
178  case kDJEpatches: if(patch->IsEMCal() || !patch->IsJetHighRecalc()) selected = false; break;
179  };
180  if(selected) result.emplace_back(patch);
181  }
182  return result;
183 }
184 
185 std::vector<const AliEMCALTriggerPatchInfo *> AliAnalysisTaskEmcalClusterNoiseTriggers::SelectFiredPatchesByTrigger(const TClonesArray &list, ETriggerThreshold_t trigger) const {
186  std::vector<const AliEMCALTriggerPatchInfo *> result;
187  EPatchType_t patchtype;
188  switch(trigger) {
189  case kThresholdEG1: patchtype = kEGApatches; break;
190  case kThresholdEG2: patchtype = kEGApatches; break;
191  case kThresholdDG1: patchtype = kDGApatches; break;
192  case kThresholdDG2: patchtype = kDGApatches; break;
193  case kThresholdEJ1: patchtype = kEJEpatches; break;
194  case kThresholdEJ2: patchtype = kEJEpatches; break;
195  case kThresholdDJ1: patchtype = kDJEpatches; break;
196  case kThresholdDJ2: patchtype = kDJEpatches; break;
197  default: return result; // unsupported patch type - return empty list
198  };
199  for(auto patch : SelectAllPatchesByType(list, patchtype)) {
200  if(patch->GetADCAmp() > fOnlineThresholds[trigger]) result.emplace_back(patch);
201  }
202  return result;
203 }
204 
205 std::vector<std::string> AliAnalysisTaskEmcalClusterNoiseTriggers::GetAcceptedTriggerClusters(const char *triggerstring) const {
206  auto clusters = PWG::EMCAL::Triggerinfo::DecodeTriggerString(triggerstring);
207  std::vector<std::string> selclusters;
208  selclusters.push_back("ANY");
209  bool isCENT(false), isCENTNOTRD(false);
210  for(const auto &c : clusters){
211  if((c.Triggercluster() == "CENT") && !isCENT) { // don't count clusters multiple times
212  selclusters.push_back("CENT");
213  isCENT = true;
214  } else if((c.Triggercluster() == "CENTNOTRD") && !isCENTNOTRD) { // don't count clusters multiple times
215  selclusters.push_back("CENTNOTRD");
216  isCENTNOTRD = true;
217  }
218  }
219  if(isCENT && isCENTNOTRD) selclusters.push_back("BOTH");
220  else {
221  if(isCENT) selclusters.push_back("ONLYCENT");
222  if(isCENTNOTRD) selclusters.push_back("ONLYCENTNOTRD");
223  }
224  return selclusters;
225 }
226 
227 int AliAnalysisTaskEmcalClusterNoiseTriggers::GetNumberNonOverlappingPatchAreas(const std::vector<const AliEMCALTriggerPatchInfo *> &firedpatches) const {
228  std::vector<const AliEMCALTriggerPatchInfo *> patchareas;
229  for(const auto patch : firedpatches) {
230  if(!patchareas.size()) patchareas.push_back(patch); // first patch always considered as new area
231  else {
232  bool overlapFound = false;
233  for(const auto refpatch : patchareas) {
234  if(!refpatch) {
235  AliErrorStream() << "Ref patch null" << std::endl;
236  AliErrorStream() << "Patchareas has size " << patchareas.size() << std::endl;
237  AliErrorStream() << "Firedpatches has size " << firedpatches.size() << std::endl;
238  }
239  if(!patch){
240  AliErrorStream() << "Test patch null" << std::endl;
241  AliErrorStream() << "Patchareas has size " << patchareas.size() << std::endl;
242  AliErrorStream() << "Firedpatches has size " << firedpatches.size() << std::endl;
243  }
244  if(HasOverlap(*refpatch, *patch)) {
245  // patch has overlap with allready accepted patch - discard
246  overlapFound = true;
247  break;
248  }
249  }
250  if(!overlapFound) patchareas.emplace_back(patch); // New non-overlapping patch found
251  }
252  }
253  return patchareas.size();
254 }
255 
256 bool AliAnalysisTaskEmcalClusterNoiseTriggers::HasOverlap(const AliEMCALTriggerPatchInfo &ref, const AliEMCALTriggerPatchInfo &test) const {
257  int testcolmin = test.GetColStart(), testcolmax = test.GetColStart()+test.GetPatchSize()-1,
258  testrowmin = test.GetRowStart(), testrowmax = test.GetRowStart()+test.GetPatchSize()-1,
259  refcolmin = ref.GetColStart(), refcolmax = ref.GetColStart()+ref.GetPatchSize()-1,
260  refrowmin = ref.GetRowStart(), refrowmax = ref.GetRowStart()+ref.GetPatchSize()-1;
261  if((InRange(testcolmin, refcolmin, refcolmax) && InRange(testrowmin, refrowmin, refrowmax)) ||
262  (InRange(testcolmax, refcolmin, refcolmax) && InRange(testrowmax, refrowmin, refrowmax))) return true;
263  return false;
264 }
265 
267  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
268  if(!mgr){
269  std::cerr << "No analysis manager defined" << std::endl;
270  return nullptr;
271  }
272 
273  bool isAOD = (mgr->GetInputEventHandler()->IsA() == AliAODInputHandler::Class());
274 
275  std::stringstream taskname;
276  taskname << "EmcalRecalcPatches_" << suffix;
278  mgr->AddTask(task);
279  task->SetEnableSumw2(true);
280 
289 
290  // setting cluster container
292  clustercont->SetClusTimeCut(-20e-9, 15e-9); // Set time cuts to [-20,+15] ns
293  clustercont->SetExoticCut(true);
294 
295  std::stringstream outfilename, outlistname;
296  outfilename << mgr->GetCommonFileName() << ":" << taskname.str();
297  outlistname << "Histos_" << taskname.str();
298  mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
299  mgr->ConnectOutput(task, 1, mgr->CreateContainer(outlistname.str().data(), TList::Class(), AliAnalysisManager::kOutputContainer, outfilename.str().data()));
300 
301  return task;
302 }
std::vector< TString > fSelectedTriggers
! Triggers selected for given event
std::vector< std::string > GetAcceptedTriggerClusters(const char *triggerstring) const
std::vector< const AliEMCALTriggerPatchInfo * > SelectAllPatchesByType(const TClonesArray &patches, EPatchType_t patchtype) const
static AliAnalysisTaskEmcalClusterNoiseTriggers * AddTaskEmcalClusterNoiseTriggers(const char *suffix)
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.
AliClusterContainer * AddClusterContainer(const char *n)
Create new cluster container and attach it to the task.
void FillTHnSparse(const char *name, const double *x, double weight=1., Option_t *opt="")
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
bool HasOverlap(const AliEMCALTriggerPatchInfo &ref, const AliEMCALTriggerPatchInfo &test) const
int GetNumberNonOverlappingPatchAreas(const std::vector< const AliEMCALTriggerPatchInfo * > &diredpatches) const
Analysis of high- tracks in triggered events.
TClonesArray * fTriggerPatchInfo
!trigger patch info array
std::vector< const AliEMCALTriggerPatchInfo * > SelectFiredPatchesByTrigger(const TClonesArray &patches, ETriggerThreshold_t trigger) const
void test(int runnumber=195345)
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.
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.
void SetClusTimeCut(Double_t min, Double_t max)
static TString ClusterContainerNameFactory(Bool_t isAOD)
Get name of the default cluster container.