AliPhysics  8d00e07 (8d00e07)
AliAnalysisTaskEmcalJetConstituentQA.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 <cmath>
30 #include <iomanip>
31 #include <iostream>
32 #include <sstream>
33 #include <string>
34 #include <THistManager.h>
35 #include <TLinearBinning.h>
36 #include <TLorentzVector.h>
37 #include <TVector3.h>
38 
39 #include "AliAnalysisManager.h"
41 #include "AliAODInputHandler.h"
42 #include "AliClusterContainer.h"
45 #include "AliEmcalJet.h"
46 #include "AliInputEventHandler.h"
47 #include "AliJetContainer.h"
48 #include "AliLog.h"
49 #include "AliTrackContainer.h"
50 #include "AliVCluster.h"
51 #include "AliVEvent.h"
52 #include "AliVTrack.h"
53 
57 
58 using namespace EmcalTriggerJets;
59 
62  fHistos(nullptr),
63  fNameTrackContainer(""),
64  fNameClusterContainer(""),
65  fTriggerSelectionString(""),
66  fUseTriggerSelection(kFALSE),
67  fNameTriggerDecisionContainer("EmcalTriggerDecision")
68 {
69  this->SetUseAliAnaUtils(true);
70 }
71 
73  AliAnalysisTaskEmcalJet(name, true),
78  fUseTriggerSelection(kFALSE),
79  fNameTriggerDecisionContainer("EmcalTriggerDecision")
80 {
81  this->SetUseAliAnaUtils(true);
82 }
83 
85  if(fHistos) delete fHistos;
86 }
87 
90 
91  TLinearBinning binningz(50, 0., 1), multbinning(51, -0.5, 50.5), binningnef(50, 0., 1.), binningR(50, 0., 0.5), binningptconst(200, 0., 200.), binningptjet(20, 0., 200.),
92  binningNCell(101, -0.5, 100.5), binningFracCellLeading(100, 0., 1.), binningM02(100, 0., 1.);
93 
94  const TBinning *jetbinning[4] = {&binningptjet, &binningnef, &multbinning, &multbinning},
95  *chargedbinning[7] = {&binningptjet, &binningnef, &multbinning, &multbinning, &binningptconst, &binningz, &binningR},
96  *neutralbinning[8] = {&binningptjet, &binningnef, &multbinning, &multbinning, &binningptconst, &binningptconst, &binningz, &binningR},
97  *binningHighZClusters[7] = {&binningptjet, &binningnef, &binningptconst, &binningz, &binningNCell, &binningFracCellLeading, &binningM02};
98 
99  fHistos = new THistManager(Form("histos_%s", GetName()));
100  for(auto c : fNamesJetContainers){
101  auto contname = dynamic_cast<TObjString *>(c);
102  if(!contname) continue;
103  fHistos->CreateTHnSparse(Form("hJetCounter%s", contname->String().Data()), Form("jet counter for jets %s", contname->String().Data()), 4, jetbinning);
104  fHistos->CreateTHnSparse(Form("hChargedConstituents%s", contname->String().Data()), Form("charged constituents in jets %s", contname->String().Data()), 7, chargedbinning);
105  fHistos->CreateTHnSparse(Form("hNeutralConstituents%s", contname->String().Data()), Form("neutral constituents in jets %s", contname->String().Data()), 8, neutralbinning);
106  fHistos->CreateTHnSparse(Form("hHighZClusters%s", contname->String().Data()), "Properties of high-z clusters", 7, binningHighZClusters);
107  }
108 
109  for(auto h : *(fHistos->GetListOfHistograms())) fOutput->Add(h);
110  PostData(1, fOutput);
111 }
112 
115  if(!tracks) tracks = GetParticleContainer(fNameTrackContainer);
116  const auto clusters = GetClusterContainer(fNameClusterContainer);
117  if(fNameTrackContainer.Length() && !tracks){
118  AliErrorStream() << "Track container " << fNameTrackContainer << " required but missing ..." << std::endl;
119  return kFALSE;
120  }
121  if(fNameClusterContainer.Length() && !clusters){
122  AliErrorStream() << "Cluster container " << fNameClusterContainer << " required but missing ..." << std::endl;
123  return kFALSE;
124  }
125 
126  // Event selection
127  AliDebugStream(1) << "Trigger selection string: " << fTriggerSelectionString << ", fired trigger classes: " << fInputEvent->GetFiredTriggerClasses() << std::endl;
128  if(fTriggerSelectionString.Contains("INT7")){
129  // INT7 trigger
130  if(!(fInputHandler->IsEventSelected() & AliVEvent::kINT7)) return false;
131  } else if(fTriggerSelectionString.Contains("EJ")){
132  auto triggerclass = fTriggerSelectionString(fTriggerSelectionString.Index("EJ"),3); // cleanup trigger string from possible tags
133  AliDebugStream(1) << "Inspecting trigger class " << triggerclass << std::endl;
134  // EMCAL JET trigger
135  if(!fMCEvent){ // Running on Data
136  if(!(fInputHandler->IsEventSelected() & AliVEvent::kEMCEJE)) return false;
137  if(!fInputEvent->GetFiredTriggerClasses().Contains(triggerclass)) return false;
138  }
139 
141  auto triggerdecisions = dynamic_cast<PWG::EMCAL::AliEmcalTriggerDecisionContainer *>(fInputEvent->FindListObject(fNameTriggerDecisionContainer.Data()));
142  if(!triggerdecisions) {
143  AliErrorStream() << "No offline trigger selection available" << std::endl;
144  return false;
145  }
146  else if(!triggerdecisions->IsEventSelected(triggerclass.Data())) return false;
147  }
148  } else return false;
149 
150  AliDebugStream(1) << "Event is selected" << std::endl;
151 
152  for(auto jc : fNamesJetContainers){
153  auto contname = dynamic_cast<TObjString *>(jc);
154  if(!contname) {
155  AliErrorStream() << "Non-string object in the list of jet container names" << std::endl;
156  continue;
157  }
158  const auto jetcont = GetJetContainer(contname->String().Data());
159  if(!jetcont){
160  AliErrorStream() << "Jet container with name " << contname->String() << " not found in the list of jet containers" << std::endl;
161  continue;
162  }
163  AliDebugStream(2) << "Reading " << jetcont->GetArray()->GetName() << std::endl;
164 
165  for(auto jet : jetcont->accepted()){
166  AliDebugStream(3) << "Next accepted jet, found " << jet->GetNumberOfTracks() << " tracks and " << jet->GetNumberOfClusters() << " clusters." << std::endl;
167  Double_t pointjet[4] = {std::abs(jet->Pt()), jet->NEF(), static_cast<double>(jet->GetNumberOfTracks()), static_cast<double>(jet->GetNumberOfClusters())},
168  pointcharged[7] = {std::abs(jet->Pt()), jet->NEF(), static_cast<double>(jet->GetNumberOfTracks()), static_cast<double>(jet->GetNumberOfClusters()), -1., 1., -1.},
169  pointneutral[8] = {std::abs(jet->Pt()), jet->NEF(), static_cast<double>(jet->GetNumberOfTracks()), static_cast<double>(jet->GetNumberOfClusters()), -1., 1., -1., -1.},
170  pointHighZCluster[7] = {std::abs(jet->Pt()), jet->NEF(), -1., -1., -1., -1., -1.};
171  fHistos->FillTHnSparse(Form("hJetCounter%s", contname->String().Data()), pointjet);
172  TVector3 jetvec{jet->Px(), jet->Py(), jet->Pz()};
173  if(tracks){
174  for(decltype(jet->GetNumberOfTracks()) itrk = 0; itrk < jet->GetNumberOfTracks(); itrk++){
175  const auto trk = jet->TrackAt(itrk, tracks->GetArray());
176  if(!trk) continue;
177  if(trk->Charge()){
178  pointcharged[4] = std::abs(trk->Pt());
179  pointcharged[5] = std::abs(jet->GetZ(trk));
180  pointcharged[6] = jet->DeltaR(trk);
181  fHistos->FillTHnSparse(Form("hChargedConstituents%s", contname->String().Data()), pointcharged);
182  } else {
183  // particle level jets
184  pointneutral[4] = pointneutral[5] = std::abs(trk->E());
185  pointneutral[6] = std::abs(jet->GetZ(trk));
186  pointneutral[7] = jet->DeltaR(trk);
187  fHistos->FillTHnSparse(Form("hNeutralConstituents%s", contname->String().Data()), pointneutral);
188  }
189  }
190  }
191  if(clusters){
192  for(decltype(jet->GetNumberOfClusters()) icl = 0; icl < jet->GetNumberOfClusters(); icl++){
193  const auto clust = jet->ClusterAt(icl, clusters->GetArray());
194  if(!clust) continue;
195  TLorentzVector ptvec;
196  clust->GetMomentum(ptvec, this->fVertex, AliVCluster::kHadCorr);
197  pointneutral[4] = std::abs(clust->GetHadCorrEnergy());
198  pointneutral[5] = std::abs(clust->GetNonLinCorrEnergy());
199  pointneutral[6] = jet->GetZ(ptvec.Px(), ptvec.Py(), ptvec.Pz());
200  pointneutral[7] = jetvec.DeltaR(ptvec.Vect());
201  fHistos->FillTHnSparse(Form("hNeutralConstituents%s", contname->String().Data()), pointneutral);
202 
203  if(pointneutral[6] > 0.95) {
204  pointHighZCluster[2] = pointneutral[4];
205  pointHighZCluster[3] = pointneutral[6];
206  pointHighZCluster[4] = clust->GetNCells();
207  pointHighZCluster[5] = *std::max_element(clust->GetCellsAmplitudeFraction(), clust->GetCellsAmplitudeFraction()+clust->GetNCells());
208  pointHighZCluster[6] = clust->GetM02();
209  fHistos->FillTHnSparse(Form("hHighZClusters%s", contname->String().Data()), pointHighZCluster);
210  }
211  }
212  }
213  }
214  }
215 
216  return kTRUE;
217 }
218 
220  using AnalysisHelpers = EMCalTriggerPtAnalysis::AliEmcalAnalysisFactory;
221  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
222  if(!mgr) {
223  std::cerr << "[AliAnalysisTaskJetConstituentQA::AddTaskEmcalJetConstituentQA(EE)] No analysis manager provided ..." << std::endl;
224  return nullptr;
225  }
226 
227  std::stringstream taskname;
228  taskname << "constituentQA_" << trigger;
229  auto task = new AliAnalysisTaskEmcalJetConstituentQA(taskname.str().data());
230  task->SetTriggerSelection(trigger);
231  mgr->AddTask(task);
232 
233  auto inputhandler = mgr->GetInputEventHandler();
234  auto isAOD = false;
235  if(inputhandler->IsA() == AliAODInputHandler::Class()) isAOD = true;
236 
237  TString tracksname, clustername;
238  AliParticleContainer *tracks(nullptr);
239  AliClusterContainer *clusters(nullptr);
240  if(!partmode) {
241  tracksname = AnalysisHelpers::TrackContainerNameFactory(isAOD);
242  tracks = task->AddTrackContainer(tracksname);
243  task->SetNameTrackContainer(tracksname);
244  tracks->SetMinPt(0.15);
245 
246  clustername = AnalysisHelpers::ClusterContainerNameFactory(isAOD);
247  clusters = task->AddClusterContainer(clustername);
248  task->SetNameClusterContainer(clustername);
249  clusters->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
250  clusters->SetClusHadCorrEnergyCut(0.3);
251  } else {
252  tracksname = "mcparticles";
253  tracks = task->AddParticleContainer(tracksname);
254  task->SetNameTrackContainer(tracksname);
255  tracks->SetMinPt(0.);
256  }
257 
258  // create jet containers for R02 and R04 jets
259  std::array<double, 4> jetradii = {{0.2, 0.3, 0.4, 0.5}};
260  for(auto r : jetradii) {
261  std::stringstream contname;
262  contname << "fulljets_R" << std::setw(2) << std::setfill('0') << int(r*10.);
263  auto jcont = task->AddJetContainer(AliJetContainer::kFullJet, AliJetContainer::antikt_algorithm, AliJetContainer::E_scheme, r, AliJetContainer::kEMCALfid, tracks, clusters, "Jet");
264  jcont->SetName(contname.str().data());
265  task->AddNameJetContainer(contname.str().data());
266  jcont->SetMinPt(20.);
267  }
268 
269  std::stringstream contname, outfilename;
270  contname << "JetConstituentQA_" << trigger;
271  outfilename << mgr->GetCommonFileName() << ":JetConstituentQA_" << trigger;
272  if(partmode) {
273  contname << "_part";
274  outfilename << "_part";
275  }
276  mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
277  mgr->ConnectOutput(task, 1, mgr->CreateContainer(contname.str().data(), AliEmcalList::Class(), AliAnalysisManager::kOutputContainer, outfilename.str().data()));
278 
279  return task;
280 }
EMCal fiducial acceptance (each eta, phi edge narrowed by jet R)
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.
double Double_t
Definition: External.C:58
Class creating a linear binning, used in the histogram manager.
AliJetContainer * GetJetContainer(Int_t i=0) const
TObjArray fNamesJetContainers
Names of the connected jet container.
void SetUseAliAnaUtils(Bool_t b, Bool_t bRejPilup=kTRUE)
TCanvas * c
Definition: TestFitELoss.C:172
Interface for binnings used by the histogram handler.
Definition: TBinning.h:21
Container for particles within the EMCAL framework.
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
void FillTHnSparse(const char *name, const double *x, double weight=1., Option_t *opt="")
static AliAnalysisTaskEmcalJetConstituentQA * AddTaskEmcalJetConstituentQA(const char *trigger, bool parmode=kFALSE)
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
Collection of helper functions used to configure the analysis.
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
AliEmcalList * fOutput
!output list
Double_t fVertex[3]
!event vertex
AliTrackContainer * GetTrackContainer(Int_t i=0) const
TString fNameTriggerDecisionContainer
Name of the trigger decision container.
Base task in the EMCAL jet framework.
Container class for histograms.
Definition: THistManager.h:99
void UserCreateOutputObjects()
Main initialization function on the worker.
void SetDefaultClusterEnergy(Int_t d)
Bool_t fUseTriggerSelection
Use trigger selection in addition to trigger string.
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.
Container structure for EMCAL clusters.
void SetClusHadCorrEnergyCut(Double_t cut)